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Preface 


The  purpose  of  this  study  was  to  survey  and  use  variance 
reduction  techniques.  To  achieve  this,  extensive  literature 
search  was  done.  For  each  technique,  the  basic  idea,  analyt¬ 
ical  formulation,  method  of  implementation,  fields  of  applica¬ 
tion,  advantages,  disadvantages  and  any  other  specific  charac¬ 
teristics  were  identified,  presented  and  clarified.  To 
illustrate  each  technique,  all  the  above  characteristics 
were  tabulated.  Numerical  examples  with  computer  programs 
were  given  and,  finally,  examples  of  application  of  the  most 
commonly  used  techniques  were  presented. 

This  work  should  result  in  better  understanding  of 
variance  reduction  techniques  so  that  one  can  use  them  more 
efficiently. 

I  would  like  to  thank  my  thesis  advisor.  Professor 
Albert  H.  Moore,  and  my  reader.  Major  Joseph  W.  Coleman,  for 
their  continuous  patience  and  assistance.  I  also  wish  to 
thank  Ms.  Sharon  Gabriel  for  her  excellent  typing  of  this 
thesis. 


Mohamed  Refat  Elhefny 
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Abstract 


The  objective  of  this  study  was  to  find  out  about, 
analyze  and  illustrate  the  characteristics  of  variance 
reduction  techniques.  Extensive  literature  review  was 
done  to  identify  the  basic  idea,  theoretical  foundation, 
procedure  for  implementation,  fields  of  application,  and 
other  specific  characteristics  of  each  technique. 

Examples  were  given  to  show  how  to  implement  the 
commonly  used  techniques.  Computer  programs  were  written  to 
perform  those  examples.  Results  were  used  to  compare  the 
efficiency  of  different  techniques.  Three  studies  in  the 
fields  of  inventory,  queuing  and  computer  performance 
measurement  were  discussed  where  different  variance  reduction 
technqiues  were  employed.  Conclusions  and  recommendations 
were  given. 


VARIANCE  REDUCTION  TECHNIQUES  WITH  APPLICATIONS 


I. 


Introduction 


v\V 


•  v  • 


During  the  early  days  of  simulation  (1940-1960) ,  when 
computer  speeds  were  much  slower,  investigators  found 
themselves  in  a  position  where  it  was  very  expensive  to 
decrease  the  variation  of  estimates  by  increasing  the  sample 
size.  Consequently,  interest  grew  in  developing  sample¬ 
estimating  procedures  that  could  either  increase  the  pre¬ 
cision  of  estimates  for  a  fixed  sample  size  or,  conversely, 
decrease  the  sample  size  required  to  obtain  a  fixed  degree 
of  precision.  Those  procedures  are  often  referred  to  as 
Variance  Reduction  Techniques  (VRTs) .  The  underlying 
principle  in  those  procedures  is  the  utilization  of  knowledge 
about  the  structure  of  the  model  and  properties  of  the  input 
data  to  change  or  distort  the  original  problem  so  that 
special  techniques  can  be  used  to  obtain  the  desired  esti¬ 
mates  at  a  lower  cost. 

Historically,  most  of  the  underlying  statistical 
approaches  used  in  VRTs  had  been  in  use  much  earlier  for 
different  purposes,  but  during  the  period  of  1940-1960 
the  techniques  had  been  refined  for  specific  use  as  vari¬ 
ance  reduction  techniques  in  computer  simulation.  As  the 
computer  speeds  increased,  the  interest  in  those  techniques 
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declined.  However,  the  recent  increase  in  complexity  of 
computer  simulation,  due  to  handling  complex  models  for 
large  problems,  renewed  the  interest  in  the  use  and  develop 
ment  of  variance  reduction  techniques.  In  some  models  of 
complex  systems,  obtaining  a  single  sample  may  require  a 
great  deal  of  computer  time  when  a  high  speed  computer  is 
used.  In  such  cases,  the  use  of  variance  reduction  tech¬ 
niques  is  vital. 

Kahn  was  one  of  the  first  pioneers  who  clarified  most 
of  the  techniques.  He  explained  and  illustrated  most  of 
them  in  the  report  published  by  Rand  in  1956  (19) .  He 
presented  several  examples  pertaining  to  the  area  of 
radiation  transport  to  demonstrate  the  applicability  of 
VRTs.  Hammer sley  and  Handscomb  presented  the  general 
Monte  Carlo  concepts  and  methods  (11) .  The  most  compre¬ 
hensive  overview  of  the  use  of  VRTs  is  presented  in  their 
book  and  also  in  the  book  by  Spanier  and  Gelbard  (30) , 
where  standard  variance  reduction  techniques,  along  with 
several  applications  to  radiation  transport  problems,  are 
discussed.  Other  books  (16;  28)  give  less  rigorous 
summaries  of  VRTs  which  are  helpful  for  the  understanding 
of  the  basic  ideas  behind  each  technique.  On  the  other 
hand,  many  articles  (7;  22;  23),  reports  (6;  8;  12;  13), 
and  studies  (5;  10)  have  been  devoted  to  development  and 
application  o  xn  VRTs  for  a  specific  kind  of  problem. 


To  compare  different  variance  reduction  techniques 
when  they  are  applicable  to  a  certain  problem,  efficiency 
of  simulation  in  estimating  parameters  is  used.  It  is 
first  suggested  by  Hammersley  and  Handscomb  (11)  and  is 
defined  as 


Efficiency 


1 


variance  X  work 


They  also  defined  the  relative  efficiency  of  simulation 
when  applying  two  Monte  Carlo  techniques  as  the  ratio  of 
their  efficiencies. 

This  implies  that  a  reduction  in  variance  of  estimator 
is  not  worthwhile  if  the  work  required  to  achieve  it  is 
excessive.  Therefore,  one  should  take  into  consideration 
the  cost  or  the  work  required  to  achieve  the  anticipated 
variance  reduction.  In  reality,  one  cannot  estimate  the 
required  work  or  the  potential  variance  reduction  for  a 
given  method.  The  analyst  can  use  his  experience  and 
intuition  to  choose  the  suitable  method  to  solve  his 
problem.  In  some  cases,  the  use  of  any  of  the  techniques 
is  infeasible  or  improf itable,  but  if  applicable  and 
properly  used,  VRTs  can  provide  a  tremendous  increase  in 
the  efficiency  of  the  simulation. 

Shannon  (26)  stated  that  variance  reduction  techniques 
are  not  new,  but  they  are  not  widely  practiced  in  spite  of 
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the  great  saving  of  work  or  reduction  of  variance  which 
can  be  achieved  when  suitable  VRTs  are  applied  to  certain 
cases. 

The  lack  of  attention  given  to  those  techniques  is 
a  consequence  of  the  shortage  of  text  books  discussing 
them  and  the  inconvenience  to  analysts  when  using  them. 

This  effort  is  devoted  to  treat  those  problems  with  the 
hope  of  making  VRTs  more  convenient  to  use. 

In  the  next  chapter,  each  of  the  known  variance 
reduction  techniques  is  illustrated.  The  types  of 
problems  which  can  be  handled  by  the  techniques  are 
discussed. 

In  the  third  chapter,  selection  and  implementation 
of  VRTs  are  discussed  in  detail,  and  procedure  for  selection 
of  the  suitable  VRT  to  a  certain  type  of  problem  is 
illustrated.  A  lengthy  table  of  the  characteristics  of 
all  the  available  standards  is  given  where  the  description, 
criteria  of  application,  advantages,  disadvantages  and 
fields  of  application  of  each  of  the  techniques  are 
condensed  to  help  in  the  selection  of  the  suitable  tech¬ 
nique. 

In  Chapter  Four,  steps  for  implementation  of  each 
technique  are  given  in  a  simple  form  including  formulae 
for  calculating  the  estimator  and  the  variance  of  the 
estimation.  Simple  examples  for  most  of  the  standard 


techniques  are  used  to  demonstrate  the  implementation  of 
the  techniques. 

In  the  fifth  chapter,  selected  examples  of  real 
applications  of  various  VRTs  are  given.  Applications  in 
the  fields  of  inventory  control  simulation  models  (6) , 
queuing  simulation  models  (10) ,  and  computer  performance 
measurements  (7)  are  demonstrated  where  the  most  commonly 
used  VRTs  are  applied. 


Variance  Reduction  Techniques 


Classification  of  Variance  Reduction  Techniques 


All  VRTs  are  concerned  with  increasing  the  accuracy 
of  Monte  Carlo  estimates  of  parameters  at  a  fixed  sample 
size  or  decreasing  the  sample  size  required  to  achieve  a 
certain  degree  of  accuracy.  In  general,  VRTs  aim  to  improve 
the  efficiency  of  the  simulation  process  when  contrasted 
with  crude  (direct  or  straightforward)  Monte  Carlo  simula¬ 
tion  which  attempts  to  create  true-to-life  or  actual 
modeling  for  the  underlying  process.  In  crude  Monte  Carlo 
simulation,  random  sampling,  flows  through  the  model  and 
sampling  probability  distributions  are  chosen  to  reflect 
the  real  situation  as  exactly  as  possible.  On  the  other 
hand,  VRTs  attempt  to  increase  the  effectiveness  of  Monte 
Carlo  simulation  by  one  of  the  following  approaches: 

1.  Modifying  the  sampling  process 


2.  Utilization  of  approximate  or  analytical  infor¬ 


mation 


3.  Studying  the  system  within  a  different  context. 
According  to  these  approaches,  the  known  VRTs  can  be 
classified  though  many  of  them  are  closely  related,  which 
makes  it  difficult  to  completely  classify  them.  The 
suggested  classification  of  most  of  the  known  VRTs  is 
given  as  follows  (24) : 


Modification  of  the  Sampling  Process 


Importance  Sampling 
Russian  Roulette  and  Splitting 
Systematic  Sampling 
Stratified  Sampling 

Use  of  Analytical  Equivalence 

Expected  Values 
Statistical  Estimation 
Correlated  Sampling 
History  Reanalysis 
Control  Variates 
Antithetic  Variates 
Regression 

Specialized  Techniques 

Sequential  Sampling 
Adjoint  Formulation 
Transformation 
Orthonormal  Functions 
Conditional  Monte  Carlo 

Modifying  the  sampling  process  is  usually  achieved 
by  using  more  effective  sampling  techniques  or  altering 
the  sampling  distributions.  This  approach  is  beneficial, 
if  not  necessary,  to  handle  simulation  problems  involving 
very  low  probability  events.  In  such  a  case,  a  modified 
sampling  scheme  is  required  to  increase  the  number  of 
occurrences  of  these  rare  events. 

Using  the  analytical  equivalence  is  another  approach 
for  reducing  the  variance  of  estimation  in  the  simulation 
process.  Since  analytical  procedures,  if  available,  are 
usually  preferable  to  simulation,  one  should  replace  the 
results  obtained  through  simulation  at  any  part  of  the 


process  by  the  available  analytical  results  or  estimates. 

In  addition  to  sampling  modification  and  analytical 
equivalence,  there  are  certain  specialized  techniques  that 
can  be  used  to  achieve  variance  reduction.  These  procedures 
may  include  the  application  of  one  or  more  of  the  above 
techniques . 

Presentation  of  Variance  Reduction  Techniques 

Modification  of  the  Sampling  Process.  Techniques  under 
this  class  have  several  common  characteristics  in  that  they 
all  reduce  the  variance  of  estimate  by  sampling  from  a 
probability  distribution  different  from  the  true  physical 
one.  This  will  help  by  observing  events  of  interest  more 
often  and  hence  decreasing  the  computing  time  and  effort. 
These  techniques  also  preserve  the  actual  process  of  system 
in  the  simulation  model,  while  only  the  probability  associ¬ 
ated  with  each  event  is  altered. 

Importance  Sampling 

Concepts  of  the  Technique.  In  this  method, 
sampling  is  forced  to  concentrate  in  the  more  important 
regions.  In  other  words,  probabilities  of  occurrence  of 
events  are  biased  in  a  known  fashion  so  that  the  resulting 
bias  can  be  adjusted  when  interpreting  the  results. 

The  idea  can  be  simply  illustrated  by  considering 


tossing  a  pair  of  dice.  If  one  is  interested  in  the 
occurrence  of  three  as  a  sum  of  the  two  top  faces,  one 


could  bias  each  die  toward  the  numbers  one  and  two  in  a 


known  fashion.  The  computation  of  the  results  should  be 
altered  according  to  the  information  from  the  biasing 
scheme  to  unbias  the  answers. 

Mathematically,  the  importance  sampling  can  be 
illustrated  by  considering  the  Monte  Carlo  estimate  of 
parameter  I  where 

I  *  E  [g  (x)  ]  =  /  g(x)  f(x)  dx 

The  crude  Monte  Carlo  procedure  for  estimating  I  would 
be  as  follows: 

1.  Select  a  random  sample  xir  ...,xN  from  the 
distribution  with  density  function  f (x) 

2.  Estimate  I  using 

l  N 

1  =  N  J1 

The  sample  variance  of  this  estimate  is  given  by 


Considering  another  distribution  f*(x)  ,  one  can  write 

I  -  /  g-(-j  -f  W  f *  (x)  dx  (2.4) 


<  w 


where  f*(x)  t  0  when  g(x)f(x)  t  0  . 


If  we  sample 


from  f*(x)  taking  x^..^  randomly,  a  new  estimator 
can  be  calculated  as 


/v 


1  N 
±  V 

N  i=l 


g<x.)f  (x.) 
f*(xi) 


(2. 


f  (x±) 

Each  sample  should,  then,  be  weighed  by  in  the 

f*(x.) 

final  result.  This  variance  reduction  procedure  will  lead 
to  a  sample  variance  given  by 


.  N  g  (x . )  f  (x . ) 

•i  *  i£r  <sr  *  *!> 


f*(xi) 


(2 


Consider  the  expected  value  EC (1^  -  I) 2 □  that  is 


g  (xi)f  (x±) 
f*(xi) 


i)  2D 


(2 


If  f*(x)  =  ?  M  an<a[  f*(x)  is  non-negative,  then 
EC (Ij  -  I)2D  =  0  which  is  a  desirable  situation  that  can 
be  achieved  only  if  I  is  known  and  f * (x)  is  chosen  to 
be  equal  to  flx^  .  Since  I  is  always  unknown,  one 

has  to  use  available  information  about  the  problem  to 
choose  f *  (x)  as  close  to  ?  (2F.1  as  possible.  If  one 

fails  to  choose  a  suitable  f*(x)  to  sample  from, 
importance  sampling  can  give  a  worse  result  than  that  of 


crude  Monte  Carlo;  that  is,  when 


ECs2  -  S*D  =  J  g2  (x)Cl  -  f  (x)dx 

1  f*(x) 


is  not  positive. 


Importance  Sampling  for  More  than  One 
Variable.  If  the  functions  f  and  g  are  functions  of 
a  vector  of  random  variables  £  ,  one  can  take  a  random 

sample  •••  from  a  selected  probability 

function  f*(X)  and  estimate  the  parameter  I  as 


1  «  g(^A)f 


-  ij 


i=l  f(Xi) 


with  the  variance 


'l  ”  N-l  i£1  ^ 


.  N  g  (X . )  f  (X  . ) 
A  V  r  A  1 


|2  - 


f*(xi) 


(2.8) 


(2.9) 


Also  in  this  case,  if  we  select  f  *  (x)  *  % llllSlZl  ,  then 
ECdj^  -  I)21  will  be  zero  indicating  the  best  selection 
of  f * (x)  which  cannot  be  done  without  knowing  I  ,  the 
parameter  to  be  estimated.  Due  to  multidimensionality  of 
the  function  f(x)  ,  it  is  difficult  to  develop  f*(x) 
efficiently.  Therefore,  a  conditional  importance  function 
may  be  selected  instead  of  f*(2{).  In  the  simplest  case. 
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when  only  two  variables  x,y  are  in  the  function  f,g  , 
the  parameter  I  can  be  estimated  as 


I  *  /  g(x,y) f (x,y)dxdy  =  /  g (x,y) f (x) f (y/x) dxdy 

x,y  x,y 


J  g(x,£)f  (x)  f*(x)f  (y/x)  dxdy 
x,y  f* (x) 


(2.1 


So  one  can  select  X^,X2»...,XN  randomly  from  a  function 

f* (x)  and  Y, ,Y_, . . . ,Y„  from  f (y/X . ) f * (X . )  and  estimate 
X  i  N  XX 


I  as 


;  1  .  9<xi'vi>f(Xi> 

1  ‘  "  t*  (X±) 


The  sample  variance  will  be 


(2.1 


„  g(X.,Y.)f  (X.) 

A  — — — -n2  - 1?) 

N  1  f*(xi)  1 


(2.1 


In  this  case 


eCi,  -  I)2:  =  {/  /  g2  (x,y) f  (y/x) dydx  -  I2} 

A  N  x  f*(x)  y 


Using  the  relation 


ECg2 (x,y)/x3  =  /  g2 (x,y)f (y/x)dy 

y 

one  can  specify  the  best  importance  function  to  be 


f*(x)  =  f  <x>  {ECg2 (x,y) /x}*8 


/{E[g2  (x^J/xDJ^f  (x)dx 


which  will  reduce  EC (1^  -  I)23  to  be 


EC  (1^  -  I)23  =  ^  {/{ECg2  (x,y)/x)  }**  f(x)dx32} 


If  we  have  more  than  two  variables,  the  best  importance 
function  can  be  expressed  as 

f  (x){E[g2  (x^J/x]}1* 


/{ E[g2  (x,y)/xDPf  (x) dx 


The  vector  y  stands  for  all  the  random  variables  except 
x  .  The  estimate  of  I  and  the  variance  of  sample  are 
expressed  in  a  way  similar  to  that  of  the  two  dimensional 
case.  It  should  be  noticed  that  the  selection  of  the  best 
f * (x)  can  be  done  only  if  we  know  the  estimator  I  for 
which  the  whole  simulation  is  made,  so  one  can  only  select 
a  good  f*(x)  guided  by  the  given  formula  for  best  f * (x) 


Russian  Roulette  and  Splitting 

Concepts  of  the  Technique.  Von  Neumann  and  Ulan 
first  used  these  two  techniques  in  particle  diffusion 
problems.  The  combination  of  the  two  techniques  can  be 
extremely  effective  when  some  knowledge  about  the  importance 
of  the  regions  of  the  distribution  is  available.  If  the 
problem  is  structured  as  a  series  of  events  that  can  be 
examined  at  various  stages,  at  some  of  these  stages,  one 
can  tell  whether  a  process  would  contribute  to  the  desired 
result  or  not.  If  the  state  of  a  certain  stage  is  not  of 
interest,  the  process  will  be  killed  off  with  a  known 
probability.  This  is  called  Russian  Roulette.  On  the 
other  hand,  if  the  process  is  an  interesting  state,  addi¬ 
tional  investigation  might  be  conducted  by  increasing  the 
number  of  simulations  starting  from  that  situation.  This 
is  called  Splitting.  As  mentioned,  the  combined  technique 
can  be  very  effective  in  multistage  problems  such  as  random 
walk,  subsystems  in  series,  etc.  It  could  also  be  useful 
in  simulations  involving  a  large  number  of  discrete  situa¬ 
tions  such  as  queuing  systems  in  which  large  numbers  of 
individuals  are  being  tracked.  In  such  systems,  Russian 
Roulette  can  be  used  to  decrease  the  number  of  individuals 
being  tracked  by  removing  an  individual  at  a  certain  stage 


in  the  problem  with  probability  .  Otherwise,  that 

individual  is  allowed  to  continue  in  the  system  with  a  new 


weight  W  =  (1  -  P^)  =  1/q^  (24) .  This  can  be  repeated 

for  the  other  individuals  and  so  the  number  of  individuals 
in  the  system  can  be  decreased.  Splitting  can  increase  the 
number  of  tracked  individuals  in  the  system  by  replacing 
an  individual  with  weight  V?  by  n  individuals;  each  of 
them  has  the  weight  W/n  .  Those  individuals  can  indepen¬ 
dently  proceed  through  the  system  keeping  their  assigned 
weights. 


Application  of  the  Technique  to  a  Two- 
Stage  Problem.  In  a  two-stage  proble,  if  X  is  the  random 
observation  from  the  first  stage  and  Y  is  that  from  the 
second,  the  estimation  of  I  =  E[g(x,y)D  can  be  calculated 


*  =  j19(xi'Yi) 


(2.18) 


where  a  sample  of  pairs  of  values,  (X^Y^,  (X2,Y2)  ,  ...  , 
(XN,  Y^)  are  generated  from  the  given  distribution  of 
x  and  y  .  If  some  values  of  X  would  lead  to  more 
interesting  results  than  others,  one  can  use  Russian 
Roulette  and  Splitting  to  divide  the  states  in  the  first 
stage  into  the  two  following  sets: 


S^  :  The  set  of  states  which  will  be  terminated 
by  Russian  Roulette  with  probability  P  =  1-q  ,  but  if  the 


simulation  is  continued  for  a  state,  the  estimated  parameter 
will  by  weighted  by  1/q 


S2  :  The  set  of  states  which  will  continue  in 
the  simulation.  Each  will  be  split  into  n  simulations 

w 

with  weight  =  —  for  each, 

n 

The  procedure  would  be  repeated  for  N  starting 
simulations  and  the  modified  estimator  will  be 


I  = 


g(xi,Yi) 


xi£Si 


VS2 


n  g(X.  ,Y  ) 


(2.19) 


which  is  the  unbiased  estimator  for  I  (24) . 

The  sample  variance  in  this  case  is  given  by 


N 


s 2  _  N  {-7  I2  -  I2} 

H-1  *N  ii1  1  1 


(2.20) 


n 

where  I.  *  0  ,  g(X.,Y.)/q  or  £  g(X.,Y.)  according 

i  ii  j=1  i  i 

to  contribution  to  the  estimator  from  the  history  i  ,  and 


I 


1  N 
-  y 

N  * 


i=l 


(2.21) 


Weight  Standard  for  General  Application  of 
the  Technique.  If  the  problem  to  be  simulated  is  broken 
into  N  regions,  two  weights  W„  and  W  will  be 

H  »  JL  • 


ii-n 


assigned  to  each  region  i  .  When  a  history  enters  region 


i  ,  its  current  weight  W  will  be  compared  to  the  region 
weight  standard  using  the  following  rules  (24) : 


as  follows: 


If  W  <  WL  ,  Russian  Roulette  is  applied 


*  kill  the  history  with  P  =  1  - 


*  the  history  will  survive  with 
W 

P  =  sj —  and  its  new  weight  will  be  W. 


2.  If  W  >  WR  ,  Splitting  is  applied  as  follows 
i 

*  find  n  such  that  W  -  n  W_  <  W_ 

*i  Ai 

*  create  n  histories  which  starts  from 
this  point,  each  with  weight  W 


W  -  n  WT 


*  with  probability 


,  create 


one  more  history  starting  from  the  same 
point  with  weight  WA 


3.  If  W_  <  W  <  W„  ,  let  the  history  continue 
Li  Hi 

in  the  simulation  without  any  change. 

The  above  procedure  will  be  used  under  the  assumption 

of  approximately  constant  importance  for  each  region.  The 

importance  of  a  region  is  inversely  proportional  to  its 

average  weight  W.  .  This  means  that  histories  moving 

Ai 
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into  a  region  of  higher  importance  (lower  weight)  will  be 
split,  while  those  moving  to  a  region  of  lower  importance 
(higher  weight)  will  suffer  Russian  Roulette.  To  increase 
the  efficiency  of  computer  time  utilization,  a  fixed  weight 
should  be  used  for  all  histories  in  a  region  of  constant 
importance.  The  high  and  low  weight  standards,  Wu  ,  WT  , 
are  used  only  to  define  the  upper  and  lower  limits  for 
triggering  Russian  Roulette  and  Splitting  processes.  They 
should  be  used  only  when  another  VRT  is  used  besides  Russian 
Roulette  and  Splitting  (24) . 

Selection  of  the  Suitable  Criteria.  There 

are  three  parameters  from  which  one  should  choose:  weight 

standards,  probability  of  kill  and  number  of  splitted 

histories.  The  best  selection  is  the  one  which  minimizes 

the  variance  in  estimate.  It  is  difficult  to  perform  this 

optimum  selection  so  the  results  from  importance  sampling 

analysis  can  help  where  the  weight  standards  for  a  given 

—  k 

region  will  be  proportional  to  (ECg2(x)3)  ^  which  means 
that  the  weight  standards  should  be  high  in  regions  of  low 
value  and  low  in  regions  of  high  value. 

Systematic  Sampling 

Concept  of  the  Technique.  Systematic 


sampling  is  a  structured  modification  of  sampling  procedure 
to  reduce  the  variance  of  estimation  of  the  parameter.  This 


technique  always  results  in  variance  reduction  without 
involving  any  significant  additional  effort,  so  it  should 
be  applied  when  that  is  possible,  in  spite  of  the  small 
improvement  it  provides.  The  technique  is  applicable  to 
any  Monte  Carlo  problem  which  has  a  probability  distribu¬ 
tion  to  characterize  the  initial  conditions.  There  are 
two  methods  to  implement  systematic  sampling  in  Monte 

oo 

Carlo  technique  to  estimate  the  parameter  I  =  /  g(x)f(x)dx 

—  00 

with  a  reduced  variance. 

Method  I .  The  range  of  the  density  function 
f (x)  is  divided  into  N  regions  with  equal  areas;  each 
equals  where  N  is  chosen  between  5  and  50.  It  is 
clear  that 

i  =  /  f  (x)  dx  ,  i=l,  ...,  N  (2.22) 

XeLi 

where  is  the  length  of  the  ith  interval  (region) .  If 

a  sample  of  random  numbers  R^,  ...  ,  R  is  selected  from 
U(0,1)  ,  the  following  sequence  of  numbers  will  be 

generated: 

R . 

Rij  =  TT  +  '  i=1'2'--*'n  '  Hm.mN 

For  each  value  of  i  ,  this  procedure  assigns  a  value 
R. .  to  each  interval  j  ,  then  a  corresponding  value  of 


the  random  variable  is  determined  from 


-  /  f(x) dx  f  i  1 #  • • • i n  g  3”^# • • • |N 

_  nn 


The  parameter  I  can  then  be  estimated  as 


-Mi- 

n  i=i  1 


(2.2 


where 


1  N 

[i  =  .1  g(X±j)  ,  i  =  1# * •  »n 


(2.2 


The  sample  variance  in  this  case  is 


S2  =  -2-  {-l  l2.  -  I2} 

n-1  n  l 


(2.2 


Notice  that  in  this  case  we  generate  only  n  random  numbers 
from  U(0,1)  and  then  generate  n  X  N  realizations  on  the 
range  of  f  (x)  i'“- 


Method  II.  This  method  is  generally  better 
than  the  first  one  to  perform  systematic  sampling.  In 
this  case  n  independent  samples  are  allocated  to  each  of 
the  regions.  This  is  done  by  selecting  R^j  ,  i=l,...,n  , 
j=l f . . . >N  from  U(0,1)  ,  then  n  random  numbers  are 

allocated  to  each  of  the  N  regions  using  the  relation 
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3  -  R. 


Rij 


N  3  '  1  1"**'n  ' 


and  the  corresponding  realization  X^j  is  determined  from 

X.  . 

.  13 

R±i  -  /  f(x)dx  (2. 


The  estimator  I  and  the  variance  of  the  sample,  S2  , 
are  determined  using  the  same  formulae,  Eqs.  2.23  and  2.25. 

The  second  method  will  always  give  better  variance 
reduction,  although  it  requires  larger  number  of  random 
samples  from  U(0,1)  .  In  the  two  cases,  the  variance 

reduction  is  approximately  proportional  to  N2 


Stratified  Sampling 

Stratified  sampling  is  similar  to  systematic 
sampling,  but  better  efficiency  is  achieved  by  taking  more 
samples  from  the  region  of  larger  variance.  It  is  a  way 
to  combine  the  features  of  systematic  sampling  with  those 
of  importance  sampling.  It  can  be  considered  as  a  special 
case  of  systematic  sampling  where  the  optimum  distribution 
of  samples  among  the  regions  is  attempted.  Usually, 
systematic  sampling  and  stratified  sampling  can  handle  the 
same  type  of  problems,  but  the  latter  is  recommended  when 
additional  information  is  available  about  region  contribu¬ 
tions  to  the  total  variance.  In  that  case,  additional 
reduction  in  the  variance  can  be  achieved. 
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Structure  of  Stratified  Sampling  Scheme. 
The  range  of  f (x)  is  broken  up  into  N  regions  of 
length  1^,  L2,...,LN  ,  respectively.  The  length  of  the 

jth  region  is  selected  in  accordance  with  a  specified 
probability 


P •  =  J  f (x)  dx  ,  j=l,...,N  (2. 

J  XeLj 

Notice  that,  if  Pj  =  ,  j=l,...,N  ,  the  same  sampling 

structure  of  systematic  sampling  will  be  obtained  where 
N  regions  of  equal  areas  are  used.  The  rule  to  select 
Pj  is  to  select  such  that  the  variance  in  g(x)f(x)  is 
the  same  in  each  interval.  After  determining  the  lengths 
of  the  intervals,  the  numbers  of  samples  from  each  interval 
nj  ,  j— 1, . . . ,N  should  be  determined.  If  the  total  number 
of  samples  is  n  where 


N 

n  *  In. 
j=l  3 


An  unbiased  (24)  estimate  for  I  is 


I 


N 

I 

j=l 


P  . 

_1 

nj 


n . 

C  t 

i=l 


(2. 


N 

I  p .  I . 
.  jil  3  D 


(2. 
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where 


^  I  9<V 


"J  i-1 


(2.31) 


The  sample  variance  in  this  case  can  be  estimated  as 


S2  = 


N  P.2 

J. 


xu 


l  -Jt  l  Cg(x.  .) 
jii  V1  i-i  v 


1 i * 


N  n.P. 


n  a  •  r  .  1  j  a  9 


(2.32) 


As  in  the  case  of  systematic  sampling,  the  stratified 

sampling,  when  compared  with  crude  Monte  Carlo,  has  an 

2 

efficiency  proportional  to  N 

Selection  of  the  optimum  number  of  samples  from  each 
internal  n^  is  a  difficult  task.  Consider 


ec  (i-i)2:  = 


N 

e[(  i  p.i.-i)2: 
1-1  3  11 


N  P  .  20 . 2 

l  -J— i_ 
j-1  n  j 


(2.33) 


where  o.2  is  the  variance  in  the  j  interval 

l 
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rav 


(2.3 


=  n.  ECd.  -  l.)^ 


If  n^'s  are  selected  to  minimize  (2.33)  subjected  to 
(2.34),  the  n^'s  should  be  selected  to  satisfy 


n  P.  o . 


y  p.  o. 

j-i 3  3 


(2.3 


Notice  that  o/  s  are  not  known,  but  they  can  be  estimated 


using 


«  l  Cg(xij)  -  i  y 

n>l  i=l  13  3 


=  Jj_  c  -L  I  g*  (X. .)  -  1*3 
nt-1  nj  i=i  13  3 


(2.3 


where  n^  samples  are  arbitrarily  selected  in  each  interval 
An  iterative  scheme  can  be  structured  to  estimate  n^ 

Analytical  Equivalence  Technique 


VRTs  in  this  group  are  based  on  using  prior  knowledge 
of  the  processes  involved  to  form  analytical  or  approximate 
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solutions  to  the  problem  being  simulated.  This  can  also 
mean  that,  if  one  can  find  a  related  process  which  can  be 
exactly  solved  using  analytical  or  other  low  variance 
techniques,  he  can  derive  the  difference  between  the  exact 
and  related  processes  using  Monte  Carlo  technique.  Many 
of  the  techniques  under  this  group  are  very  closely 
related  in  the  principles  and  ideas  involved. 

Use  of  Expected  Value.  This  method  is  based  on 
the  fact  that  an  analytical  determination  of  parameter 
estimator  is  usually  preferred  to  the  results  of  simulation 
procedures.  Since  Monte  Carlo  estimation  of  a  parameter 
is  an  estimation  of  its  expected  value,  the  technique  is 
so  called  and  it  is  applicable  where  the  expected  value  of 
portions  of  the  model  can  be  determined  analytically  with¬ 
out  losing  an  essential  element  of  the  simulation.  Expected 
value  method  can  be  used  in  multistage  problems  where  the 
expected  value  of  the  parameter (s)  can  be  analytically 
determined  in  one  stage  or  more.  For  example,  consider  the 
two  stage  problem  where  X  is  selected  from  f (x)  at 
the  first  stage  and  Y  is  selected  from  f(y/x)  at  the 
second.  Repeating  the  process  N  times,  crude  Monte 
Carlo  estimation  of  I  is 

~  .  N 

I  -  £  l  g(X.,Y  •)  (2.37) 

N  i=l  1  1 
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If  it  is  possible  to  determine  analytically  f(y/x)  in  the 
second  stage  for  a  given  value  X  from  the  first  stage, 
the  simulation  process  can  be  simplified  where  N  samples 
from  f  (x)  ,  are  generated  and  the  estimator  for 

I  can  be  calculated  as 

*  ,  N 

!e  =  £  E[g(y/Xi)]  (2. 

which  is  an  unbiased  estimator  of  I  since  ECl„D  =  I 

£ 

The  sample  variance  in  this  case  is  given  by 

s’  “  X  e’Cg.y/x,  -  i’3  «• 

In  addition  to  the  simplification  of  the  simulation 
processes,  the  above  technique  always  gives  better  results 
compared  with  crude  Monte  Carlo. 

It  should  be  noticed  that  it  is  not  always  possible 
to  calculate  the  expected  value  of  a  portion  of  the 
simulation  analytically.  An  approximation  of  the  expected 
value  obtained  by  another  variance  reduction  technique  may 
be  used.  In  some  cases,  a  portion  of  simulation  cannot  be 
replaced  by  its  expected  value  even  if  it  is  analytically 
determined.  In  those  cases,  the  second  and  higher  moments 
may  be  important  in  the  simulation  procedures  and  not  only 
the  expected  value. 
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Statistical  Estimation.  In  this  technique, 


the  stochastic  process  is  not  removed  from  the  simulation, 
but  the  expected  value,  rather  than  the  simulation  result, 
is  used  in  the  estimation. 

If  one  step  in  a  simulation  is  a  random  choice  between 
reaching  some  final  outcome  or  continuing  in  the  simulation 
process,  then  Statistical  Estimation  can  be  used.  In  crude 
Monte  Carlo,  a  random  number  R  would  be  generated  at  this 
step  and  if  R  <  P(Y^/J)  ,  then  the  history  would  be 

terminated  with  score  1.  If  R  >  P(Y^/X)  ,  then  the 

history  would  continue  with  no  score  being  made.  After 
N  histories,  the  estimate  for  probability  of  reaching 
Yf  would  be 


W  -  & 


(2.40) 


where  n  is  the  number  of  histories  terminated  at  Y^ 

In  statistical  estimation,  the  same  simulation  process 
is  used,  but  the  estimation  technique  is  changed.  Every 
time  the  particular  step  is  encountered,  a  contribution  of 
P(Y^/X)  is  added  to  the  estimate,  regardless  of  the  actual 
outcome  of  the  simulation.  The  final  estimate  is  then 
given  by 

PSE<Vf>  -  s  Xl  P<VXij>  (2-41> 
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where  the  second  summation  runs  over  all  occurrences  of 
the  possibly  final  step  in  the  course  of  jth  simulation. 
An  estimate  of  the  variance  may  be  calculated  from 


where 


=  n  [i  ?  ;2 . 

5PI  lN  . 1 


'“se'V3 


(2.42) 


Pi  -  l 


(2.43) 


The  use  of  statistical  estimation  will  always  improve  the 
variance  of  estimation,  but  it  can  be  particularly  useful 
if  the  probability  of  reaching  the  desired  end  point  is 
small  in  all  intermediate  stages.  It  becomes  essential 
when  the  probability  of  the  end  point  becomes  vanishingly 
small.  If  there  were  many  intermediate  stages  which 
could,  with  very  low  probability,  reach  the  desired  end 
point,  then  statistical  estimation  might  calculate  the 
desired  result  with  good  accuracy. 


\m 


Correlated  Sampling 

Concept  of  the  Technique.  Correlated 
sampling  can  be  one  of  the  most  powerful  VRTs  due  to  the 
wide  applicability  of  the  technique,  as  well  as  the  large 
efficiency  gains  which  can  be  obtained.  If  the  primary 
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objective  of  a  simulation  is  to  determine  the  effect  of 
a  small  change  in  the  system,  crude  Monte  Carlo  approach 
would  make  two  independent  runs,  with  and  without  the 
change  in  the  system,  then  subtract  the  results  obtained. 
Usually,  the  difference  will  be  smaller  than  either  of  the 
two  outputs,  but  the  variance  of  the  difference  will  be 
the  sum  of  the  variances  in  the  two  runs.  In  such  cases, 
the  use  of  correlated  sampling  can  be  essential  to  obtain 
statistically  significant  results.  If  the  two  simulations 
use  a  common  random  number  at  comparable  stages  in  the 
computation,  the  correlation  in  results  in  the  case  of 
correlated  sampling  will  reduce  the  variance  of  estimation 
much  more.  Another  way  of  viewing  correlated  sampling 
through  random  number  control  is  to  realize  that  the  use 
of  the  same  random  numbers  will  generate  identical  histories 
in  those  parts  of  the  system  which  are  the  same,  so  that 
the  difference  in  results  will  be  due  to  the  difference 
in  the  two  systems.  This  will  increase  the  efficiency 
compared  to  uncorrelated  cases.  Correlated  sampling  can 
be  utilized  in  the  following  types  of  simulations: 

*  calculation  of  the  effect  of  small  change 
in  the  system 

*  difference  in  parameter  in  two  or  more 
similar  cases  is  of  more  interest  than  its  absolute  values 


in  them 


* 


performance  of  parametric  study  of  several 


similar  problems 

*  the  answers  to  unknown  problems  is  to  be 
estimated  using  the  known  answer  of  a  similar  problem. 

Analytical  Formulation  of  Correlated 
Sampling.  Let  the  integrals  1^  ,  Ij  characterize  two 
different  but  related  problems. 

1 1  =  I  f  (x)  g^xjdx 


and 


r2  =  /  f2(y)  s2(y)dy 

If  the  main  interest  is  the  difference 


crude  Monte  Carlo  approach  will  perform  two  separate 
simulations  where,  in  the  first,  the  estimator  of  1^  j 
calculated  as 


A 


N 


1  " 


using  a  sample  X^,...,X  selected  randomly  from  f^(X) 
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and 

I2  is  estimated 

as 

A 

*2  = 

1  N 

N  ,|,  g2(Yi) 

using  a  sample  Y,,.. 

.,Yn  selected  randomly  from  f 2 (y) 

The 

difference  is  then  estimated  as 

/V  A  A 

A  = 

The 

variance  in  this 

case  is 

A 

o2 (A)  =  a 

*  (I,)  +  a2  (I2)  -  2  CovU,,^) 

where 

=  EC (1,-1,) 2: 

2 

a2(I2) 

=  eC(i2-i2)2: 

and 

fM 

<  H 

*» 

•H 

<  M 

> 

0 

u 

A 

=  e[ (1,-1,) (i2-i2)] 

=  eC  (i,,i2) D  -  i,i2 

A  A 

Now  if  I,  ,  I2 

are  positively  correlated,  then 
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0 


(2.54) 


Cov 


‘W 


and  the  variance  in  the  correlated  case  will  be  less  than 
that  realized  with  no  correlation. 

History  Reanalysis 

Concept  of  the  Technique.  History  reanalysis 
is  essentially  a  form  of  correlated  sampling,  except  that 
one  does  not  actually  run  a  second  simulation  using  the 
same  random  numbers  as  in  the  first.  Instead,  the  results 
of  the  first  simulation  are  reanalyzed  to  calculate  the 
answer  for  the  second  process  (24) .  This  technique  reduces 
the  variance  due  to  correlation  and  cuts  down  the  computa¬ 
tional  time  involved  since  the  second  simulation  is  not 
actually  performed.  The  technique  can  handle  the  same 
types  of  problems  listed  in  the  correlated  sampling  case 
with  the  condition  that  the  differences  in  the  systems 
being  simulated  must  be  expressible  as  a  difference  in 
probability  distribution  or  in  the  scoring  function. 

Analytical  Formulation.  Assume  that  there 
are  two  problems  which  involve  estimating  1^,  I2  as 
given  by  Eqs.  (2.44)  and  (2.45).  Assume  also  that  a  random 
sample  X^,...,XN  has  been  obtained  from  f^(X)  .  The 

estimator  for  1^  is  as  usual 


11-27 


*1  =  N  J,  *l<Xi> 


(2.55) 


Since 


I  =  /  f9(X)  g_(x)dx  = 


g, (x) f _ (x) 

/  Vx)  f ,  cx>  dx 


an  estimator  for  I2  can  be  obtained  using 


(2.56) 


i  “  WW 

2  '  N  iii  £i(Xi> 


(2.57) 


where  f^(X^)  ^  0  whenever  g2  (X^)  f  2  (X^)  ^  0  .  The 

sample  variance  of  I2  is 


N  ,1  ?  r92«i>f  2(Xi> 

H-l  ‘N  ^  E7TXT1 


(2.58) 


To  calculate  the  effect  of  correlation,  it  is  necessary  to 
estimate  the  variance  of  the  difference  directly.  That  is, 
if 


Ai  " 


g2  (xi)  f2  (xi) 

W 


-  *i(xi> 


(2.59) 


is  the  difference  in  the  ith  history  and 
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A 


1  N 
—  j 
N  i=l 


(2.60) 


is  the  average  difference,  then  the  sample  variance  is 


.2  _ 


N 


N 


n  ri  r  u  :2-, 
i= T  CN  2  Ai  ‘  A  11 


i=l 


(2.61) 


Control  Variates 

Concept  of  the  Technique.  In  many  simulation 
problems,  there  exist  simplifications  or  approximations  to 
the  problem  having  analytic  or  closed  form  solutions.  In 
these  cases,  the  analytic  information  can  be  used  to  reduce 
variance  by  what  is  referred  to  as  control  variates.  In 
this  technique,  the  difference  between  the  problem  of 
interest  and  some  analytical  models  approximating  it  is 
simulated.  The  gain  in  variance  reduction  or  estimating 
accuracy  is  proportional  to  the  degree  of  correlation 
between  the  true  process  and  the  analytical  model  used. 

This  approach  has  a  wide  range  of  applicability  and  it  is 
very  useful  when  analytical  representations  of  simplified 
models  exist  (24). 

Analytical  Formulation.  Consider  the 

integral 


oo 

I  =  /  g(x)f(x)dx 

—00 


(2.62) 
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Assume  that  it  is  possible  to  get  a  function  h(x)  whose 
expected  value  is  known  or  can  be  analytically  determined, 
and  which  approximates  g(x)  .  If  the  value  of  the 
integral 


oo 

0  =  |  h  (x)  f  (x)  dx 

—00 


(2.63) 


is  known,  then  the  integral  I  can  be  expressed  as 


+  00  00 

I  =  /  h(x)f(x)dx  +  J  Cg(x)  -  h(x)]f(x)dx 

—  00  —  oo 

oo 

=  6  +  /  Cg(x)  -  h(x)D  f  (x)  dx  =  0  +  I..  (2.64) 

—  00  •*- 


The  function  h(x)  is  called  the  control  variate  for 
g(x)  .  Since  8  is  known  or  can  be  calculated  analyt¬ 
ically,  simulation  is  needed  for  estimation  of  1^  , 


00 

1 2  =  /  Cg(x)  -  h(x)D  f(x)dx 

—  00 


(2.65) 


This  can  be  performed  using  crude  Monte  Carlo  by  selecting 
a  sample  X^,...,XN  from  f(x)  and  using 


A 


.  N  .  N 

N  9<Xi>  '  N  h<Xi» 


1 

N 


N  „ 


i=l 


Ai 
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(2.66) 


vV  s'  v*V 


where 


A±  =  g(Xi)  -  h(Xi) 


An  estimate  of  the  sample  variance  is 


s2  - 


(2.67) 


(2.68) 


where 


A  = 


=  ±  l  A. 

N  i=l  1 


(2.69) 


The  efficiency  of  the  control  variates  technique  depends 
on  the  degree  of  similarity  between  g(x)  ,  h(x) 


Antithetic  Variates 

Concept  of  the  Technique.  This  technique  is 
similar  to  the  control  variates  approach,  except  that  a 
negatively  correlated  function  is  chosen.  This  negative 
correlation  is  used  to  reduce  the  variance  of  estimation. 
Another  difference  between  control  variates  and  antithetic 
variates  is  that  the  expectation  of  the  chosen  function 
need  not  be  known.  Antithetic  variates  approach  can  be 
implemented  in  several  methods.  Two  of  these  methods  are 
discussed  here. 
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Method  I .  If  the  parameter  I  where 


00 

I  -  /  g  (x)  f  (x)  dx 

00 


(2. 


is  to  be  estimated  using  antithetic  approach,  an  unbiased 

A 

estimator  1^ 


-  It  <2- 


is  found  using  crude  sampling.  A  second  unbiased  estimator 

A  A  A 

I2  for  I  is  selected  such  that  1^  ,  I ^  are  negatively 

A  A 

correlated.  Linear  combination  of  I1  ,  I2  is  a  third 
unbiased  estimator  for  I  which  can  be,  for  example 


6  =  a  11  +  (1  -  a)  I2 


(2. 


In  this  method,  «<  is  chosen  to  be  simply  i,  then 


A  A 


e  -  f  +  V 


(2. 


will  be  an  unbiased  estimator  for  I  with  a  variance 
given  by 


(6)  =  \  (1^  +  \  a2  (I2)  + 


1 

1 


(2.74) 
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Since  Ij  is  chosen  to  be  negatively  correlated  with  1^ 


then 


cov  d1,i2)  -  0 


This  could  make  the  variance  of  the  combined  estimator  6 
smaller  than  the  variance  of  either  of  the  two  estimators 

A  A. 

I1  '  I2  * 

A  convenient  way  to  accomplish  this  method  is  to 
generate  a  set  of  random  numbers  from  U(0,1) 

and  two  negatively  correlated  sets  of  random  variables 
X1#...,XN  and  X^,...,X^  are  obtained  using  the  same 
set  of  random  numbers  where#  for  each  selected  random 
number  R.^  ,  the  corresponding  X^  ,  X^  are  calculated 

from 


R.  =  J  f  (x)  dx 
—  00 


-  R^  =  /  f  (x)  dx 


The  negative  correlation  between  each  pair  of  values 

A  A 

X .  ,  X'  is  clear;  then  the  two  estimators  I.  and  I. 
11  12 

will  be  negatively  correlated.  Defining 


6i  =  ?  Cg(Xi)  +  g(X^)D 
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the  estimator  of  I  using  antithetic  variates  will  be 


N 


-  y 

N  , 1 


e. 


i=l 


1 

2N 


N 

I  Cg  (x  )  +  g  (x^ )  3 

i=l  1  1 


(2.79) 


The  sample  variance  is  given  by 


S2 (6)  = 


N 

N-l 


<ir 


N 

l 

i=l 


et  -  e2} 


(2.80) 


Method  II.  In  this  method  of  implementing 
the  antithetic  approach,  we  try  to  find  the  value  of  a 
which  best  improves  the  estimation  efficiency.  This 
approach  can  be  viewed  as  a  combination  of  antithetic 
variates  and  stratified  sampling  where  the  range  of  f (x) 
is  divided  into  two  strata,  <  x  <  X„  and  Xu  <  x  <  ® 

M  M 

If  a  random  number  Ri  is  selected  from  U(0,1)  ,  a  pair 

of  values  X.^  ,  X^  can  be  calculated  from 


and 


aR.  =  /  f(x)dx 

—00 

X' 

a+(l  -  a)R^  =  /  f(x)dx 


(2.81) 

(2.82) 


which  means  that  X^  is  selected  from  the  range 

<  x  <  XM  and  X".  from  XM  <  x  <  «  .  Also,  X.  ,  X' 
are  negatively  correlated.  In  this  case,  the  combined 
variable  will  be 
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(2.83) 


ei  =  a  g(Xi)  +  (1  -  a)  g(X^) 


and  the  new  unbiased  estimator  for  I  is 


e 


.  N  .  N 

s  I  6.  =  i  l  Cag(X  )  +  (l-a)g(x;)]  (2.84) 

N  i«l  1  N  i=l  1  1 


with  the  sample  variance 


•  2  _ 


N 


=  JL _/!  y  0- 

N-l  {N  i£1  ei 


-  0  2  } 


(2.85) 


If  a  is  properly  selected,  this  method  can  give  much 
better  results  than  the  first  simple  one.  Selection  of 
a  is  a  difficult  task,  but  a  rule  of  thumb  is  to  select 
a  such  that 


g(XM)  *  a  g(XL)  +  (l-oOgUy)  (2.86) 

where  Xy  and  XL  are  the  upper  and  lower  limits  of  the 
range  f (x)  .  Alternatively,  a  can  be  determined  using 

a  trial  and  error  method  to  obtain  the  optimum  efficiency. 

Regression.  Regression  techniques  can  be  applied 
to  a  wide  variety  of  Monte  Carlo  simulations  to  produce 
unbiased  estimators  for  a  set  of  parameters  (integrals)  when 
correlation  between  them  is  known  to  exist.  Regression 
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technique  will  make  use  of  this  correlation  to  reduce  the 
variance  of  estimation. 


Analytical  Formulation.  If  a  set  of 
integrals  I^f...,I  are  to  be  estimated,  regression 
can  be  applied  to  determine  the  minimum  variance  unbiased 

A  A 

estimators  by  a  set  of  estimates  0^,....,0  (n  Z  P)  such 
that 


E  ( 0  . )  -  a • •  I.  +  ....  +  a.  I  ,  j— l,...,n  (2.87) 

J  J  ^  ^  Jr  r 


where  a^  ,  j=l,...,n  and  i=l,...,P  is  a  set  of  known 
constant.  If  the  coefficients  a^  in  Eq  (2.87)  is 
represented  in  matrix  form. 


all 

• 

• 

,alP 

A  = 

a21 

• 

• 

a2P 

(2.88) 

_anl 

• 

* 

anP 

and  a  sample 

consisting 

of 

N  independent 

sets  of  simulated 

values  for 

' 

then  one 

calculates 

A 

1 

N 

9  . 
3 

_  ± 
N 

l 

i=l 

e . . 

13 

,  i=  1 ,  • 

• « f  n 

(2.89) 

to  construct  the  column  matrix 
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(2.90) 


Now,  an  estimate  for  the  set  of  integral  I  where 


(2.91) 


is  given  by  I  , 


A  ,-*T  -*--1  -1  ■+-T  ->-1  * 

I  =  (A  V  A)  A  V  6 


(2.92) 


where 


1 


fry: 


r o . 
•v  ■ 


ty- 


IK' 

>  V 
ifv 


,-<«s 

& 

y  *• 


h>\ 


■W.S.' 


J.'  SM' 


V  = 


V 


11 


V 


In 


V 


21 


V 


2n 


•  • 


V 


nl 


V 


nn 


(2.93) 


^  ^  m 

is  the  covariance  matrix  for  6.,..., 6  and  A  is  the 

l  n 

transposition  of  A  .  That  is. 


Vij  =  E  C{ei  "  E(ei)}  {ej  "  E(8j)>3 


i  1 , ■ . . , n  ,  3=1 , . . . ,n 


which  could  be  estimated  as 


(2.94) 


N 


vij 


J1  (eki  -  ei»  (eM  -  V 


kj  "j 

1-1, . . . ,n  ,  3  l,...,n 


(2.95) 


where  6^  is  calculated  from  Eq  (2.90),  and 
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k>'.- 


i  Vi viV«  1  'll 


V 


11 


V 


21 


V. 


N1 


V 


In 


V 


2n 


V 


nn 


The  new  unbiased  estimator  is 


I*  =  (£T  V-1  A)-1  AT  V  e 


(2.96) 


(2.97) 


Generally,  it  is  difficult  to  formulate  the  estimators 

/v  /\ 

0l'’**'6n  *  This  limits  the  applicability  of  the  method 

in  real  situations. 

Specialized  Techniques 

This  group  of  VRTs  includes  those  techniques 
which  are  useful  to  a  specific  kins  of  problem.  Some  of 
those  techniques  are  not  well  developed  or  extremely 
specialized.  However,  some  of  those  techniques  are  the 
only  way  to  get  a  considerable  variance  reduction  in 
certain  cases.  In  this  section,  the  most  common  tech¬ 
niques  of  the  above  characteristics  will  be  discussed. 
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XT 


Conditional  Monte  Carlo.  In  some  cases. 


is  is  hard  to  estimate  the  parameter 


e  =  E  E 4>  (a)  □ 


(2.98) 


where  a  is  a  random  vector  distributed  over  a  space 
A  with  a  probability  density  function  f(a)  .  This  may 
be  due  to  the  complexity  of  f(a)  .  One  way  to  deal 
with  such  a  problem  is  to  embed  the  space  A  in  a 
product  space  C  =  A  x  B  where  B  is  suitably  chosen. 
Each  point  in  C  can  be  expressed  as 


c  =  (a,b) 


and  a  can  be  considered  as  a  function  of  c  which 
maps  the  points  of  C  to  A  .  If  we  sample  a  random 
vector  y  =  (a,  8)  from  C  with  a  probability  density 
function  h(c)  ,  a  mapping  of  y  to  a  is  obtained 


(2.99) 


which  is  a  random  vector  of  A 


In  general,  a  will 


not  have  the  desired  density  function  f  ,  so  an  appro¬ 
priate  weighting  function  should  be  used  to  compensate 
for  that. 

If  we  choose  g(c)  =  g(a,b)  ,  an  arbitrary  real 
function  defined  on  C  such  that 
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.V  /- 


V  V  *-  *.  %  \  * 


and 


G(a)  =  /  g(a,b)  db  ^  0  for  any  (a,b) 

B 


h(c)  ^  0  for  any  (c) 


a  suitable  weighting  function  could  be 


w(c)  =  f  (a)  g  (c)  J (c)  /G(a)h(c) 


where  J (c)  is  the  Jacobian  of  the  transformation 
c=(a,b)  which  can  be  written  as 


J(c)  =  J(a,b) 


dadb 

dc 


Since  a  is  the  first  coordinate  of  c  ,  the 
identity  holds  (11) 

/  <j>  (a)  f(a)da  =  /  da  $  (a)7-ffa-}  J  g(a,b)db 

A  A  G(a)  B 


=  / 

AxB 


<t>  (a)  f  (a)  g  (c) 
G  ( a )  h  ( c ) 


h (c) dadb 


/ 

C 


()>  (a)  w(c)  h  (c) 


dadb 

J(c) 


/  $  (a)w(c)  b(c)d  (c) 
C 


(2.100) 

(2.101) 

(2.102) 

(2.103) 

following 


(2.104) 
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This  shows  that,  if  a  is  the  first  coordinate  of  a 
random  vector  y  sampled  from  C  with  density  function 


h(c)  ,  then 


t  =  <}>(a)  W(y) 

is  an  unbiased  estimator  of  6  .  Here  B  ,  h  are 

chosen  arbitrarily  to  simplify  the  sampling  procedure. 

The  function  g  acts  as  an  importance  function  which 

should  be  selected  to  minimize  the  variation  in  t  , 

■  • 

and  hence  increase  the  precision  of  estimation. 

Conditional  Monte  Carlo  is  a  special  case  of  the 
above  theory  where  h(c)  is  a  given  distribution  on 
C  =  A  x  B  and  f (a)  =  f(a,bQ)  is  the  conditional 
distribution  of  h(c)  given  that  b  =  bQ  we  have 

h(c)d(c)  =  f(a,b)  ^  (b)  dadb 

where  i|>(b)  is  the  probability  density  function  of  B 
and  y  =  (a,g)  has  the  density  function  h(c)  .  In 
this  case 

J  (c)  =  h(c)  /  f  (a ,  b)  i|>(b) 

and  for  a  given  b 

o 


(2.105) 


(2.106) 


(2.107) 


J  ( a ,  bQ )  =  h(a,bQ)  /  f (a)  4(b0) 


Eliminating  f(»)  from  the  weight  function,  we  get 


(2.108) 


w(c) 


h(a,bD)J(a,b)g(a,b) 
h(a,b)  J(a|s0)^(b0)G(a) 


(2.109) 


This  leads  to  the  following  rule.  Let  y  =  (a,B) 
be  distributed  over  C  with  probability  density 
function 


Then 


h(c)  =  h(a,b) 


(2.110) 


t  =  <|)(a)w(Y) 


(2.11) 


where  w(y)  is  given  by  Eq  (2.111)  as  an  unbiased  estimator 
of  the  conditional  expression  of  41(a)  given  that  B  —  bQ  . 
It  is  clear  that  this  approach  requires  neither  sampling 
from  space  A  which  may  be  awkward,  nor  evaluation  of  the 
possibly  complicated  function  f  .  Besides,  we  can 
achieve  variance  reduction  in  the  estimation  if  g  is 
suitably  selected. 


Sequential  Sampling.  This  is  not  a 
specific  variance  technique,  but  rather  a  general  approach 
to  the  use  of  other  techniques.  It  is  useful  to  apply 
this  technique  when  there  is  little  or  no  apriori  infor¬ 
mation  about  the  expected  results  of  the  simulation.  In 
this  technique,  a  series  of  sequential  simulation  runs 
is  performed  to  reduce  the  variance  of  estimated  para¬ 
meters.  In  the  first  run,  little  or  no  reduction  variance 
is  achieved.  In  the  second  run,  parameters  estimated  in 
the  first  one  are  used  for  applying  other  VRTs  such  as 
importance  sampling,  Russina  Roulette,  splitting,  or 
stratified  sampling.  A  third  run  can  then  be  made  using 
the  improved  sampling  parameters  and  this  "self-learning" 
process  can  be  carried  out  repeatedly  with  the  efficiency 
of  sampling  improving  at  each  stage.  In  spite  of  the 
simplicity  of  this  approach,  little  work  on  sequential 
sampling  has  been  done  (2 H) .  Considering  this  technique 
a  trade-off  must  be  done  between  the  required  extensive 
computation  and  the  efficiency  gain  from  improved  sampling 
The  sequential  nature  of  this  technqiue  may  lead  to  more 
underbiased  or  overbiased  estimation,  if  the  initial 
choice  of  the  parameter  is  biased. 

Orthonormal  Functions.  This  VRT  can  be 


very  useful  when  applied  to  multidimensional  problems. 
In  this  method,  a  set  of  orthogonal  functions  over  a 


region  of  multiple  integration  is  defined,  then  a 
sampling  scheme  is  structured  to  permit  efficient 
sampling  over  this  region  from  a  joint  probability  density 
function.  The  procedures  to  implement  this  technique  are 
not  well  developed,  but  the  potential  gain  when  applying 
it  is  still  high. 

Adjoint  Method.  Frequently,  when  building 
a  simulation  model,  one  can  find  a  set  of  mathematical 
equations  which  is  adjoint  or  inverted  with  respect  to 
the  original  set.  In  such  cases,  a  solution  for  one  set 
of  equations  will  give  the  solution  for  the  second  set. 

The  basic  idea  in  applying  the  adjoint  method  as  a  VRT 
is  to  simulate  the  adjoint  set  of  equations  which  does 
not  represent  any  real  process,  but  is  easier  or  more 
efficient  to  simulate.  It  would  give  a  solution  which 
helps  in  estimating  the  original  parameter  direclty  or 
in  applying  another  VRT.  In  some  cases,  one  can  divide 
the  problem  into  two  parts;  in  one  of  them,  the  adjoint 
method  is  applied  while  direct  simulation  is  applied  in 
the  second  part. 

The  adjoint  method  has  been  exploited  very  success¬ 
fully  in  radiation  transport  problems  because  of  the 
precise  formulation  of  this  problem  as  a  linear  integral 
equation  for  which  an  adjoint  formulation  can  be  obtained 


This  technique  needs  more  investigation  and  further 


development  to  be  generally  applicable  in  simulation. 

Transformations .  This  method  is  a  special 
form  of  importance  sampling  which  is  characterized  by 
formulating  the  priori  information  about  the  process  in 
a  parametric,  closed  form  representation.  That  informa¬ 
tion  can  be  used  to  alter  the  sampling  procedure  by 
transformation.  This  method  has  been  largely  employed 
in  radiation  transport  calculations  where  the  function  of 
interest  have  an  approximately  exponential  form  (24). 
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III.  Selection  of  Variance  Reduction  Techniques 

In  the  preceding  chapter,  basic  concepts  of  different 
VRTs  were  discussed.  Those  techniques  are  not  equally 
efficient  when  applied  to  a  specific  problem.  The  selection 
of  a  promising  technique  for  a  particular  problem  can  cause 
considerable  difficulty  due  to  the  large  number  of  possibil¬ 
ities  available.  This  chapter  is  devoted  to  set  general 
rules  for  selection  of  the  appropriate  technique (s)  for  a 
certain  situation  or,  in  other  words,  to  show  where  each 
technique  can  generally  be  used.  To  achieve  this  goal,  a 
summary  of  properties  and  concepts  of  most  of  the  known 
VRTs  is  given  in  a  tabular  form  to  help  in  selecting  the 
suitable  technique  (s) . 

For  the  analyst  to  select  and  implement  an  appropriate 
variance  reduction  technique  or  techniques,  the  following 
systematic  procedure  should  be  applied: 

1.  Definition  of  the  problem  information  that  can  be 
used  as  a  basis  to  select  an  appropriate  technique (s) . 

2.  Selection  of  specific  technique (s)  that  should 
be  considered  for  a  given  problem. 

3.  Setting  of  basic  guidelines  to  implement  the 
selected  procedure. 

These  aspects  are  described  in  the  following  three 
sections. 
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Definition  of  Problem  Information 

The  efficiency  of  a  variance  reduction  technique  is 
strongly  related  to  the  efficiency  of  the  use  of  the 
known  information  about  the  problem.  Primarily,  it  is 
essential  to  characterize  the  aspects  of  the  problem  that 
can  help  to  indicate  the  fruitful  variance  reduction 
technique  for  this  problem.  Helpful  information  items 
are  organized  in  the  following  table.  These  items  are 
the  basic  items  needed  for  most  of  the  techniques. 

Table  3.1  presents  the  required  useful  information 
that  should  be  known  prior  to  the  selection  and  implemen¬ 
tation  of  the  suitable  variance  reduction  technique. 

Selection  of  Variance  Reduction  Technique (s) 

The  most  difficult  step  in  utilization  of  variance 
reduction  techniques  is  the  selection  of  the  suitable 
technique (s)  that  would  fit  the  problem  of  interest  and 
give  an  effective  variance  reduction  or  reduce  the  required 
sample  size  for  a  given  degree  of  accuracy.  This  diffi¬ 
culty  can  be  reduced  by  preparing  and  utilizing  the 
information  about  the  problem  listed  in  Table  3.1  and 
understanding  the  characteristics  of  different  available 
variance  reduction  techniques.  For  that  reason,  a  compre¬ 
hensive  summary  of  variance  reduction  techniques  is 
presented  in  this  section.  Having  the  information  about 
the  problem  under  investigation  in  mind,  and  understanding 
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2. 


0. 


TABLE  3.1 

PROBLEM  INFORMATION 

NEEDED  FOR  SELECTION  OF  SUITABLE  TECHNIQUE 


1.  Definition  of  nature  of  the  problem  relative  to: 

*  expected  values  (means,  variances,  etc.)  to  be  estimated 

*  sensitivities  or  variations  of  parameters  of  interest 

*  possible  mathematical  formulations  (integral  equations, 
expected  values,  etc.) 

*  any  sequential  characteristics,  such  as  independent  path; 
outcomes  depend  on  intermediate  step 


input  conditions  which  are  random  variables  to  be  sampled. 

Identification  of  portions  of  problem  or  parameter  to  be 

estimated  that  can  be: 

*  expressed  in  an  analytical  form  such  as  single  integral, 
multiple  integral,  differential  and/or  integral  equations 

*  solved  analytically,  such  as  expected  values,  variances, 
probabilities,  etc. 

*  represented  by  approximate ,  simplified  positively  correlated 
analytical  expressions 

*  represented  by  approximate,  simplified  negatively  correlated 
analytical  expressions 

*  established  as  relatively  not  important  to  final  outcomes 
compared  to  other  aspects  of  the  problem. 

Identification  of  variables  in  the  problem  which: 

"  are  very  important  to  the  expected  outcome 

*  are  not  expected  to  significantly  inpact  the  results 

*  are  strongly  correlated  with  other  variables. 


4.  Location  of  final  events  or  outcomes  of  the  problem  which: 

have  very  small  probabilities 
have  very  large  probabilities 

have  outcomes  relatively  insensitive  to  problem  parameters 

have  known  probabilities  of  occurrence  from  intermediate 
stages  in  the  problem 

are  linear  combinations  of  other  events  or  randan  variables 
have  known  correlation  with  other  events  or  outcomes. 


the  basic  idea  of  each  available  variance  reduction  tech¬ 
nique,  one  can  use  the  following  summary  of  the  character¬ 
istics  of  different  techniques  to  select  the  techniques 
that  could  fit  his  problem.  With  the  revision  of  the  problem 
information  and  selected  technique (s)  characteristics,  the 
most  suitable  technique (s)  could  be  identified. 

The  following  summary  of  variance  reduction  techniques' 
characteristics  includes  description,  supposed  criteria  of 
application,  advantages  and  disadvantages  of  each  technique 
(24) .  It  also  includes  the  typical  area  of  application  of 
each  technique.  It  should  be  noted  that,  in  many  cases, 
more  than  one  technique  can  be  separately  applied  to  the 
problem,  but  each  of  them  will  deal  with  the  problem  in 
different  approaches.  Also,  in  some  cases  more  than  one 
technique  may  be  applied  to  solve  one  problem  in  the  same 
time.  Each  of  these  techniques  will  be  used  in  one  stage 
of  the  problem. 

The  most  important  point  to  keep  in  mind  when  selecting 
and  implementing  one  or  more  of  the  variance  reduction 
techniques  for  a  certain  problem  is  that  the  applied  tech¬ 
nique  (s)  will  reduce  the  variance  of  only  one  parameter  or 
aspect  of  the  problem  being  simulated.  Using  variance 
reduction  techniques  designed  for  one  parameter  will 
usually  reduce  the  effectiveness  of  the  simulation  to 
estimate  other  parameters.  Therefore,  it  is  very  important 
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to  determine  all  of  the  results  which  will  be  desired  from 


the  simulation  before  searching  for  a  suitable  variance 
reduction  technique.  When  more  than  one  quantity  is  to  be 
estimated,  the  chosen  technique (s)  should  not  degrade  the 
efficiency  of  any  of  the  estimations.  In  many  situations, 
it  may  be  advisable  to  implement  a  different  variance 
reduction  technique  for  each  parameter. 
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_ Characteristics  of  Variance  Reduction  Techniques _ 

DESCRIPTION  CRITERIA 

OF  OF  TYPICAL 

TECHNIQUE  APPLICATION  ADVANTAGES  DISADVANTAGES  APPLICATIONS 


DESCRIPTION 


DESCRIPTION  CRITERIA 

OF  OF  TYPICAL 

TECHNIQUE  APPLICATION  ADVANTAGES  DISADVANTAGES  APPLICATIONS 
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queuing  problems 


DESCRIPTION  CRITERIA 

OF  OF  TYPICAL 

TECHNIQUE  APPLICATION  ADVANTAGES  DISADVANTAGES  APPLICATIONS 


DESCRIPTION 


TT".* 


IV.  Implementation  of  Variance  Reduction  Techniques 


Once  the  suitable  VRT  is  selected  to  be  applied  for 
a  specific  simulation  problem,  a  plan  for  implementing  the 
selected  technique  should  be  made.  In  many  cases,  the 
theory  given  in  the  literature  cannot  be  directly  used  to 
implement  the  selected  technique.  In  this  chapter,  general 
guidelines  to  implement  the  more  important  VRTs  are  given 
in  the  form  of  a  step-by-step  procedure.  Simple  examples, 
to  show  how  to  apply  these  steps,  are  presented  for  the 
most  commonly  used  techniques.  Computer  programs  were 
written  in  FORTRAN  V  to  accomplish  those  examples.  These 
examples  would  give  better  insight  of  the  implementation 
and  the  efficiency  of  the  techniques. 


Importance  Samplinc 


Guidelines  for  Implementation.  The  general  guidelines 
that  could  be  followed  to  implement  importance  sampling 
are  as  follow: 


1.  Express,  if  possible,  the  expected  value  being 
estimated  as 


I  =  Jg(x) f (x) dx 


where  x  is  the  random  variable  of  importance 
sampling  and  f (x)  is  its  density  function. 
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2.  Determine  the  functional  form  of  g(x)  analytically 
or  by  selecting  values  of  x  and  estimating  the 
corresponding  g(x) 

3.  Select  an  importance  function  f * (x)  which 
approximates  f{x)g(x) 

4.  Select  a  random  sample  X^,...,XN  from  f* (x) 
using  a  set  of  random  numbers  from  U(0,1) 

5.  Estimate  I  using  (2.5). 

6.  The  estimator  for  the  sample  variance  can  be 
calculated  using  (2.6). 

Example .  If  it  is  required  to  estimate  the  probability 
P(x  £  1)  when  f(x)  is  given  by 

f(x)  =  .01  exp(-.Olx)  (4.2) 

the  crude  Monte  Carlo  will  require  hundreds  of  sample 
values  from  f(x)  to  accurately  determine  P(l)  .  This 
is  due  to  the  fact  that  x  will  be  less  than  1  in  approx¬ 
imately  1/100  of  the  sample  values.  Applying  importance 
sampling  method  could  reduce  the  required  sample  size  for 
a  given  degree  of  accuracy.  Following  the  above  guidelines 
to  implement  importance  sampling,  the  value  of  P(l)  can 
be  estimated  as  follows: 

1.  Express  the  required  integral  as 


I 


(4.3) 


1 

=  /  f{x)dx 

0 

where  f (x)  is  given  by  (4.2). 

2.  Let  f *  (x)  be 


f*(x)  =  exp(-x)  (4.4) 

The  selected  importance  function  f * (x)  will 
give  the  opportunity  for  x  to  take  values  less 
than  or  equal  to  1  more  frequently. 

3.  Select  a  random  number  Ri  from  U(0,1)  ,  then 

determine  the  corresponding  Xi  using 

Xi  =  F*"1(Ri)  =  Jloge(l/(l-Ri))  (4.5) 


4.  Determine  the  values  of  f(X^)  ,  f*(X^) 

5 .  Repeat  steps  3  and  4  for  N  times  where  N  is 
the  chosen  sample  size. 

6.  Determine  the  estimator  for  P(l)  as 


P(l) 


N 


l 

i=l 


f  (X.) 

f*Tx~y 


i 


7.  Determine  the  sample  variance  using  (2.6). 


(4.6) 


The  above  steps  were  coded  in  a  computer  program  using 


FORTRAN  V.  The  code  and  the  results  are  given  in  Appendix 
A.  The  results  are  given  for  N  =  20  .  This  small  sample 

size  gave  an  accurate  estimation  of  P(l)  as  .01015  (the 
theoretical  value  is  .00995  ).  The  sample  variance  was 
.0252  .  If  crude  Monte  Carlo  were  used,  it  would  require 
hundreds  of  sample  values  to  give  the  same  accuracy. 

Russian  Roulette  and  Splitting 

Guidelines  for  Implementation.  To  apply  Russian 
Roulette  and  Splitting,  the  following  general  steps  can 
be  followed: 

1.  Determine  stages  of  the  problem  in  which  possible 
conditions  can  be  divided  into  N  regions  where 
each  of  them  contains  points  of  roughly  the  same 
importance . 

2.  For  each  region  choose  average  weight  standards, 

W_  ,  i=l  ,  N  .  This  weight  should  be  inversely 
Ai 

proportional  to  the  importance  of  the  region. 

3.  If  no  other  VRTs  are  used,  set  high  and  low  weight 

standards,  W„  ,  W_  ,  equal  to  W.  ;  other- 
Hi  Li  l 

wise,  W„  ,  W  should  be  sufficiently  spaced 
Hi  Li 

above  and  below  W  .  The  spacing  should  pre- 
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vent  any  unnecessary  Russian  Roulette  or  Splitting 


and  assign  approximately  eaual  weights  to  histories 
of  roughly  the  same  importance . 


When  a  history  arrives  at  a  particular  stage  in 
region  R^  with  weight  W  ,  carry  out  the  suit¬ 
able  action  according  to  the  following  cases: 


If  W  <  W, 


,  apply  Russian  Roulette  by 


killing  the  history  with  probability  1  - 


W 


W, 


or  letting  it  continue  in  simulation  with 


probability  W/W.  carrying  a  new  weight 


W, 


b.  If  WT  <  W  <  W„  ,  let  the  history  continue 
Li  Hi 


with  weight  W 


If  W  >  Wv 


,  carry  out  splitting  as  follows 


(i)  Determine  n  such  that 


0  <  W  -  n  W  <  WA 

i  l 


(ii)  Split  the  history  into  n  "daughter" 
histories,  starting  in  that  point  with 


weight  W 


Ai 


W-n  WA. 

(iii)  With  probability  — ^ -  ,  create  one 


W 


A. 

l 


more  daughter  history  with  weight  W 


Form  estimate  1^  for  each  history  i 


I . 


i  -  I  9<V  w( 


p,  daughter  of  i 


6.  The  final  estimate  of  I  is  given  by 


,  N  * 

-  5  I  ii 

N  i=l  1 


where  N  is  the  number  of  starting  histories. 
7.  The  variance  of  the  sample  is  given  by  (2.20). 


Example  of  Using  Russian  Roulette  in  System  Reliabilit; 
Rice  and  Moore  introduced  a  Monte  carlo  technique  for 
estimating  lower  confidence  limits  on  system  reliability 
(26) .  The  algorithm  which  they  gave  to  accomplish  crude 
Monte  Carlo  requires  a  generation  of  total  3000  random 
variates  from  the  normal  distribution  N(0,1)  .  Russian 

Roulette  can  be  applied  to  reduce  the  number  of  sample 
values  generated  from  N(0,1)  when  using  l his  algorithm. 

Case  1:  Series  Connection  of  Subsystems.  Consider 
a  system  of  three  subsystems  connected  in  series.  The 
reliability  of  each  subsystem  Pi  is  calculated  using 
pass  or  fail  test  as 


P. 

l 


1 


F. 

1 


N 


where  P^  is  the  number  of  failures  in  N  trials.  One 
can  roughly  say  that 


where  rg  is  the  system  reliability.  In  fact,  each 
is  a  random  variable  which  can  be  considered  asymptotically 
normal  (26) ,  accordingly  Rg  will  be  a  random  variable 
depending  on  P.^  ,  P2  and  P3  .  If  it  is  required  to  deter¬ 
mine  P(r  <  R  )  or  P(r  1  R  )  ,  a  modified  algorithm 

S  X  S  X 

using  Russian  Roulette  can  be  carried  out  as  follows: 

1.  Using  the  data  from  pass-fail  tests,  the  first 
estimate  of  reliability  of  each  subsystem  can 
be  calculated  as 


P 


i 


1  -  F./n. 

l  l 


( 


the  asymptotic  variance  is 


V. 

l 


where  q^  =  1  -  P^  ,  n^  being  the  number  of 
pass-fail  tests  of  the  subsystem  i 


2.  Draw  a  random  variable  from  N(0,1) 


3.  Calculate  the  second  estimate  P.  as 


4.  If 


P1  *  Rx 


,  skip  steps  5  through  8. 


5. 

Calculate 

A 

P2 

as  in  : 

steps  2 

and  3  where  i=2 

6. 

Calculate 

R12 

r12 

= 

A  A 

A  A 

P  P 

1  2 

(4.14) 

If  r12  ; 

R 

X 

,  skip 

steps 

7  and 

7. 

Calculate 

ft 

P3 

as  in 

steps 

2  and  3,  where  i=3  . 

8. 

Calculate 

rs 

as 

r 

s 

A 

/\ 

=  P 

/V 

A 

•  P 

1  2 

A 

A 

‘  P3 

(4.15) 

If  r  <  R  ,  then 

S  X 

9.  Add  1  to  the  number  of  trials  where  r  <  R 

s  x 

10.  Start  a  new  trial  at  step  2  and  repeat  for 
1000  trials. 

Notice  that  the  application  of  Russian  Roulette 
reduces  the  number  of  generated  normal  random  variates 
to  one  in  some  trials  and  two  in  other  trials.  In  some 
trials,  we  have  to  generate  three  random  variates  from 
normal  as  in  the  original  algorithm.  The  total  number  of 

AAA 

reductions  depends  on  the  values  of  Pi  '  P2  '  P3  anc*  Rx 

The  above  algorithms  have  been  coded  to  a  computer 
program  in  FORTRAN  V  which  is  attached,  along  with  the 
results,  in  Appendix  B.  The  attached  results  are  obtained 


IV-8 


when  P^  ,  an<*  P3  are  *7/  .8  and  .8,  respectively. 

and  Rx  runs  from  .3  to  .75.  The  number  of  reductions 

increases  with  the  increase  of  R  .  When  R  is  .75, 

x  x 

a  total  of  1654  out  of  3000  generations  from  normal 
N(0,1)  were  saved. 


Case  2:  Parallel  Connection  of  Subsystems.  The 
above  algorithm  can  be  used  to  handle  systems  consisting 
of  parallel  subsystems  with  some  modifications.  Suppose 
that  we  have  a  system  with  three  parallel  subsystems  which 
were  examined  using  pass-fail  tests.  If  it  is  required  to 
calculate  P(r  -  R  )  ,  one  can  use  the  same  algorithm 

given  for  the  series  case  to  apply  Russian  Roulette  to 
simulation  with  the  following  modifications: 

/v 

A 

a.  In  step  4,  we  skip  to  9  if  P.  Z.  R 

«L  X 

b.  In  step  6,  r12  is  calculated  as 


1  -  (1  -  Px) (1  -  P2) 


(4.16) 


and  we  skip  to  step  9  if  r12  -  Rx 

c.  In  step  8,  rg  will  be  calculated  as 


rs  =  1  -  (1  =  Px) (1  -  P2) (1  -  P3) 


(4.17] 


and  we  check  if  r  1  R 


d.  In  step  9,  we  add  1  to  the  number  of  trials 

where  r  1  R 
s  x 

The  computer  program  for  the  parallel  case  is  presented 

A 

in  Appendix  C,  along  with  the  results  for  equal  to 

.7,  i=l,2,3  and  Rx  runs  from  .905  to  .950  .  The  results 

showed  that  the  number  of  reductions  decreases  when  Rx 

increases.  When  Rx  was  .905,  a  total  of  573  reductions 

in  the  number  of  generated  random  variates  from  N(0,1)  out 

of  a  total  3000.  In  another  experiment,  more  reduction 

was  achieved  as  R  was  decreased. 

x 

Systematic  Sampling 

Systematic  sampling  can  be  implemented  using  the 
following  steps: 

1.  Determine  the  cumulative  function  for  f (x) 

Divide  its  range  into  N  intervals,  each  of 
width  1/N  .  N  should  be  between  5  and  50. 

2.  Generate  n  sets  of  N  random  numbers,  each 

from  U(0,1)  .  Denote  them  by  R^,...^^  ; 

R21'*“'R2N  ;  Rnl ' '  *  * '  RnN  * 

3.  Allocate  the  generated  random  numbers  into  the 

corresponding  intervals  using 
j  -  R.  . 

Rij  “  N  '  1-1 ' *  *  •  'n 

j=l,...,N  (4.18) 
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4.  Determine  X^j  corresponding  to  each  from 


X.  . 

.  *3 

Rii  =  /  f(x)dx 

•J  _ on 


(4 


5.  Estimate  the  integral  I  using  (2.23)  and  (2.24)  . 

6.  Estimate  the  sample  variance  using  (2.25). 

Example.  Suppose  it  is  required  to  estimate  the  value 

of  the  integral 


I 


(4 


using  Monte  Carlo  techniques.  In  fact,  the  value  of  this 
integral  can  be  easily  calculated,  but  it  will  be  used  to 
demonstrate  how  to  implement  several  VRTs.  Let  the 
integrand  of  (4.20)  be 


f  (x) 


(4 


The  crude  Monte  Carlo  procedure  to  estimate  I  requires 
the  generation  of  N  random  numbers  from  U(0,1)  ,  then 

determination  of  the  values  of  f(x)  at  these  points. 
That  is  because  the  range  of  integration  runs  between  0 
and  1.  If  the  number  of  generated  random  numbers  N  is 
small,  crude  Monte  Carlo  will  give  an  inaccurate  estimate 
of  I  with  a  large  variance  of  sample.  However,  one  can 


.19) 


.20) 


.21) 
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improve  the  accuracy  of  estimation  and  reduce  the  sample 
variance  by  using  one  of  the  VRTs.  Systematic  sampling 
can  be  applied  to  estimate  I  as  follows: 

1.  Divide  the  range  of  integration  (0-1)  to  N 
intervals  each  of  width  .25 

2.  Generate  n  random  numbers  from  U(0,1)  for 
each  interval,  then  allocate  them  inside  that 
interval  using  (4.18) .  Determine  the  value  of 
f(x)  at  this  point. 

3.  The  estimator  of  I  is  calculated  using  (2.23) 
and  (2.24) . 

2 

4.  The  sample  variance  S  is  calculated  using 
(2.25) . 

A  computer  program  in  FORTRAN  V  was  written  to  estimate 
the  integral  I  given  by  (4.20)  using  crude  Monte  Carlo, 
systematic  sampling  and  stratified  sampling.  The  purpose 
of  that  is  to  demonstrate  the  implementation  of  the  three 
techniques  and  to  compare  the  results  when  the  same  sample 
size  and  the  same  random  number  stream  are  used.  The 
results  for  the  three  methods  are  attached,  along  with  the 
program  code,  at  Appendix  D.  These  results  will  be  dis¬ 
cussed  in  the  next  section. 

Stratified  Sampling 

Stratified  sampling  can  be  implemented  using  the 
following  steps: 
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1.  Divide  the  range  of  the  variable  being  simulated 

into  N  intervals  of  length  L^,...,LN  .  N 

should  be  between  5  and  50.  can  be  selected 

so  that  the  variation  of  integrand  is  approximately 
the  same  inside  the  interval  j 

2.  Determine  P^  ,  the  probability  that  x  will 

be  in  as 


P.  =  /  f(x)dx  ,  j=l,...,N 

3  xeL . 

3 


3. 


Determine  the  number  of  sample  values,  n^  , 
j=l  ,  N  taken  from  each  interval  using  (2.34) 


(4.22) 


and  (2.35)  . 

4.  For  each  interval  j  ,  select  a  set  of  n^  random 
numbers  R^j  ,  i=l  ,  n^  from  U(0,1)  .  Allocate 

these  random  numbers  in  the  specific  interval  j 
and  calculate  the  corresponding  values  of 
using 


jY1  fX;Lj 

R  .  _.P_;  +  l  P  =  /  f  (x)  dx 


13  3 


P=1 


(4.23) 


5.  Estimate  I  using  (2.30)  and  (2.31). 

2 

6.  Calculate  S  using  (2.32). 


Example .  It  is  required  to  estimate  the  value  of  the 


integral  given  by  (4.20)  using  stratified  sampling.  The 


above  steps  can  be  used  to  implement  the  technique  as 
follows. 

1.  The  range  of  integration  (0-1)  was  divided  into 
N  =  4  intervals. 

2.  The  length  of  each  interval  was  chosen  to 

have  the  same  variation  of  f (x)  in  each  interval 
The  intervals  turn  out  to  be  (0-.36)  ,  (.36-. 62) 
(.62-. 83)  ,  (.83-1)  . 

3.  For  each  interval,  four  random  numbers  were 
generated  from  U(0,1)  and  allocated  inside  the 
interval.  The  value  of  f (x)  corresponding  to 
each  random  number  was  calculated. 

4.  The  estimator  of  I  inside  each  interval  j  was 
calculated  using  (2.31). 

5.  The  final  estimator  of  I  was  calculated  using 
(2.30)  . 

6.  The  variance  of  estimation  was  calculated  using 
(2.32)  . 

The  computer  program  in  Appendix  D  performs  these  steps 
to  implement  stratified  sampling.  However,  it  also  per¬ 
forms  systematic  sampling  and  crude  Monte  Carlo  as  stated 
before.  Results  for  the  three  methods  are  also  given  in 
Appendix  D.  From  these  results,  it  is  clear  that  the 
stratified  sampling  method  gives  the  best  estimation  among 
the  three  methods.  The  estimator  of  I  (=  .4153)  is*  the 
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nearest  to  the  exact  value  of  I  (=  .418)  .  Also,  the 

method  gives  the  lowest  sample  variance  (.0174)  .  In  fact, 

the  estimations  of  I  given  by  crude  Monte  Carlo,  systematic 
sampling  and  stratified  sampling  were  .3834,  .4396,  .4153  and 
the  sample  variances  were  .1198,  .0220,  .0174,  respectively. 
These  results  give  an  idea  about  the  efficiency  of  each  of 
the  three  methods,  taking  into  consideration  the  additional 
computational  effort  done  in  the  last  two  methods. 


Expected  Value 


When  expected  value  technique  is  applicable  to  a  certain 
simulation  problem,  the  procedure  to  implement  the  technique 
will  differ  according  to  the  role  of  the  random  process, 
which  can  be  replaced  by  its  expected  value,  in  the  overall 
simulation.  In  fact,  one  of  the  following  cases  will  fit 
the  problem  under  consideration.  For  each  possible  case, 
guidelines  to  use  the  technique  are  given. 

1.  If  the  process  to  be  replaced  by  its  expected 

value  is  a  selection  of  a  random  variable  y  from 
a  density  function  f (y)  ,  set 


y  =  E  [f  (y)  ] 


(4.24) 


and  continue  the  simulation. 

2.  If  the  process  represents  a  decision  between 

terminating  or  not  terminating  the  history,  let 
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the  history  continue  in  simulation  carrying  a 


.v 


jLm. 

P 


weight  W  given  by 

new 


W  =  W  . .  •  P  (4.25) 

new  old  s  '  ' 


where  Pg  is  the  probability  of  survival  of  the 
history  at  the  decision  point  and  W  .  .  and  VJ 
are  the  weights  of  the  history  before  and  after 
the  replaced  random  process.  For  any  parameter 
being  calculated,  an  estimate  for  each  history  can 
be  made  by  summing  the  contributions  from  that 
history;  that  is. 


ri  *  . 


where  is  the  weight  of  the  ith  history  at 

the  time  of  the  jth  contribution  to  the  final 
result.  The  final  estimate  for  the  parameter  is 
given  by 


(4.26) 


1  N 

=  -  y 

N  .S 


I; 


(4.27) 


i=l 


v*v 


and  the  sample  variance  is  given  by 


N 


JL.  {I  y  !2 

N-i  lN  i£1  l 


I2> 


(4.28) 
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When  the  contributions  to  a  parameter  from  a 
history  come  from  the  terminations  in  the  replaced 
process,  the  loss  of  weight  at  each  step  is  the 
proper  measure  for  the  expected  termination;  in 

A 

this  case  1^  for  the  ith  history  will  be 


A 


I  . 
1 


y  w  , ,  . .  (i-p  ) 

4  old, 13  s 


g(xij> 


(4 


where  j  denotes  the  jth  occurrence  of  the 

~  2 

replaced  event  in  the  ith  history.  I  and  S 

are  calculated  as  before  using  (4.27)  and  (4.28). 

3.  If  the  replaced  process  represents  a  decision 

between  two  branch  points,  the  history  must  be 

split  and  followed  from  this  point  as  two  separate 

histories,  one  at  each  branch  carrying  a  weight 

equal  to  the  branching  probability.  To  estimate 
2 

I  and  S  ,  formulas  identical  to  (4.26),  (4.27) 
and  (4.28)  can  be  used  where  the  contributions 
from  all  daughter  histories  are  considered. 


Statistical  Estimation 

Statistical  estimation  can  be  implemented  when 
applicable  using  the  following  steps. 

1.  Identify  the  stochastic  process  in  the  simulation 
which  has  the  desired  final  outcome  as  one  possible 


alternative . 


When  the  process  of  interest  is  encountered  in 
simulating  a  history,  score  the  following 
contribution 


g(X,Yf)  f(Yf/X) 

where  g(x,y)  is  the  function  being  integrated 
by  the  simulation,  Y^  is  the  desired  final  outcome 
of  the  process,  X  denotes  the  current  state  of 
all  other  variables  in  the  system,  and  f(Yf/X)  is 
the  conditional  probability  of  obtaining  outcome 
Yf  given  X  as  the  status  of  the  system. 

The  simulation  should  not  be  modified,  but  the 
stochastic  process  of  interest  is  modeled  by 
selecting  Y  from  (y/X) 

If  the  outcome  of  step  3  is  Yf  ,  no  additional 
scoring  is  to  be  made.  The  contribution  of  :this 
step  remains  g(x^,Y^)  f(Yf/X^) 

Estimate  the  total  contribution  of  history  i  ,  as 

*i  -  |  9<Sij'V  f  (Yf/5ci)  (4. 

where  j  runs  over  all  occurrences  of  the  particu¬ 
lar  process  being  estimated  in  the  ith  history. 
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6.  The  estimation  of  I  and  S  can  be  calculated 


using  (4.27)  and  (4.28). 

Correlated  Sampling 

If  there  are  two  similar  simulations  involving  only 
a  single  random  variable  and  it  is  desired  to  estimate 


A  = 


(4.31) 


where 


00 

1 1  ~  f  g1(x)£1(x)dx 

—CO 


(4.32) 


00 

I2  *  J  g2<y)f2(y)dy 

—00 


(4.33) 


then  correlated  sampling  can  be  implemented  in  this  case 
as  follows. 

1.  Generate  N  random  numbers  R^,...,R^  from  U(0,1)  . 

2.  Generate  a  random  sample  Xlf...,XN  from  f^tx) 
and  another  sample  of  the  same  size  Yi'Y2'**''YN 
from  f2(y)  using 


R. 


l 


/  f  (x)dx  =  /  f  2  (y)  dy 

00  —00 

i=l , . . .N 


(4.34) 


3.  Estimate  A  using 
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where 


A. 

1 


9l(Xi)  -  92,Yi> 


(4.36) 


4.  Estimate  the  sample  variance  using 


s2  .  _N_  (i  y 

b  N-l  1N  . 

1=1 


62  -  62) 


(4.37) 


Notice  that  the  variance  will  be  greatly  reduced 
if  f-L(x)  is  similar  to  f  2  (y)  and  g-^  (x)  is 
similar  to  g2  (y)  ,  that  is  because  the  two 

random  samples  X^,...,XN  and  Y^,...,YN  will 
be  highly  correlated. 

Control  Variates 

Control  variates  technique  can  be  implemented  using 
the  following  steps  to  evaluate  the  integral  I=/g (x) f (x) dx  . 

1.  Express  the  parameter  (s)  to  be  estimated  in 
integral  form. 

2.  Obtain  an  approximate  function  h(x)  for  each 
parameter  I  .  The  expected  value  of  0  of 
h(x)  should  be  known. 

3.  Estimate  the  integral  1^  , 


(4.38) 


I1  =  f  [9(x)  -  h(x)]  f  ( x )  dx 


I1  =  ^  J  l9(xi)  "  h(Xi)]  (4.39) 

where  X^,...,X^  are  a  random  sample  generated 
from  f (x) 

4.  Calculate  the  final  estimator  of  I  as 


I  =  0  +  ^  (4.40) 


(4.41) 

Example.  The  above  steps  were  used  to  estimate  the 
value  of  the  integral  given  by  (4.20).  Again,  the  same 
integral  is  used  to  compare  different  techniques.  The 
approximate  function  was  chosen  to  be 


5.  The  variance  of  estimation  can  be  calculated  as 


n^t  {ir  J1lg(xi)  "  h(xi>^  - 


N 


h  (x )  =  x 


(4.42) 


which  has  a  mean  value  . 5  at  the  range  of  integration 
(0-1)  .  The  steps  were  coded  in  a  computer  program  written 


A  A  .\  -V'  v-  v' 


■  A  .V 
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in  FORTRAN  V  which  is  attached,  with  the  results,  in 
Appendix  E.  This  method  estimated  the  integral  I  given 
by  (4.20)  as  .4255  with  sample  variance  .0016  when  the  sample 
size  was  20.  It  should  be  noted  that  the  accuracy  of  esti¬ 
mation  depends  on  how  far  the  approximate  function  h(x) 
mimics  f (x) 

Antithetic  Variates 

To  implement  this  technique,  one  should  find  two 
negatively  correlated  estimators  for  the  parameter  of 
interest.  A  linear  combination  of  these  two  estimators 
can  form  a  third  estimator  which  would  have  a  smaller 
variance  than  the  variance  of  either  of  the  original 
estimators.  Steps  to  implement  antithetic  variates 
technique  can  be  as  follows: 

1.  Put  the  parameter  to  be  estimated  in  integral 
form  as 

00 

I  =  /  g(x)  f(x)dx  (4.43) 

—  00 

2.  Select  a  value  for  the  parameter  <x  between  0,1. 

3.  Generate  a  set  of  N  random  numbers  from  U(0,1)  . 
denote  them  Ri'R2’***’RN  * 

4.  For  each  R  ,  calculate  two  negatively  correlated 

l 

random  variates  X^  ,  X^  .  This  can  be  accom¬ 
plished  using 
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a  =  /  f  (x)  dx 

—  oo 

x*: 

i 

1  -  a  R.  =  /  f  (x) dx 


5.  Calculate  the  values  of  the  new  random  variable. 


0±  =  a  g(Xi)  +  (l-a)g(X^)  ,  i=l,...,N 


6.  Estimate  I  using 


i 


7.  The  variance  of  estimation  is 


Notice  that  a  can  be  simply  taken  equal  to  1/2  or  an 
iterative  process  can  be  used  to  determine  the  optimum 
value  of  a  which  gives  the  minimum  sample  variance. 

Example .  Considering  again  the  evaluation  of  the 
integral  I  given  by  (4.20),  one  can  use  the  above  steps 
to  implement  antithetic  variates  technique  as  follows: 


1  N 
=  -  y 

N  .s 

1=1 


4.44) 

4.45) 
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2.  g(x) 


e  -  1 


3.  Let  a  =  1/2  initially 


4.  Generate  R^  from  I7  (0,1) 


5.  X.  =  R. 

l  l 


6.  x:  =  1  -  R. 

l  l 


7.  ei  =  oc  g(Xi)  +  (l-o)  g(Xj) 


8.  Repeat  steps  4  through  7  for  N  times 


9.  Estimate  I  as 


N 


-I  11  A 

6  =  4  l  e. 

N  i=i  1 


(4.50) 


(4.51) 


(4.52) 


(4.53) 


(4.54) 


The  above  steps  were  coded  in  a  computer  program 
which  is  attached  in  Appendix  F.  The  value  of  a  was 
changed  in  the  range  (.5-. 95).  The  results  are  attached, 
also  in  Appendix  F.  From  these  results  the  accuracy  of 
estimation  appeared  to  be  sensitive  to  the  change  of  a 
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Nl 
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Among  the  chosen  values  of  a  ,  the  initial  value, 
a  =  1/2  ,  gives  the  minimum  variance  equal  to  .0004  which 

indicates  that  this  method  can  achieve  a  great  gain  in  the 
efficiency  of  estimation  when  compared  with  crude  Monte 
Carlo  with  small  additional  computational  effort.  Notice 
that  minimum  variance  is  the  criteria  of  accuracy  of 
estimation  in  this  case  since  the  estimates  resulting 
from  using  antithetic  variates  approach  are  unbiased. 
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V.  Application  of  Variance  Reduction  Techniques 

The  earliest  applications  of  VRTs  were  found  in  particle 
and  radiation  transport  problems  where  very  low  probabilities 
are  involved  (20)  during  the  1940s  and  1950s.  In  those  days, 
the  objective  of  applying  VRTs  was  to  compensate  for  the  low 
speeds  of  computers  which  were  frustrating  when  dealing  with 
such  problems  (28) .  In  the  past  two  decades,  computer  speed 
has  increased  tremendously.  This  decreased  the  attention 
given  to  the  development  and  application  of  VRTs  for  a  while, 
then  an  increased  demand  for  applying  VRTs  appeared  due  to 
the  increasing  complexity  of  problem  simulations  which  would 
consume  a  great  amount  of  computer  time  or  result  in  reduced 
estimation  accuracy  if  none  of  the  VRTs  were  employed. 

Recently,  VRTs  have  found  a  wide  appxication  in  almost 
all  simulation  areas:  inventory  simulation  (6) ,  queuing  models 
simulation  (10),  network  analysis  (8),  reliability  studies, 
stationary  (13)  and  non-stat ionary  (14)  simulation  models, 
population  growth,  and  simulation  of  Markov  process  (12) . 

In  this  chapter,  examples  of  these  applications  in  the  fields 
of  inventory  simulation,  queuing  simulation  and  computer 
performance  measurement  are  presented.  The  aim  of  presenting 
these  examples  is  to  show  how  VRTs  can  be  applied  in  such 
fields. 
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Application  of  Variance  Reduction  Technique  for 
Inventory  Simulation 

Stockout  in  inventory  control  systems  should  be  rare 
events.  Therefore,  variance  reduction  is  an  important 
experimental  design  issue  to  estimate  accurately  the  mean 
values  of  those  rare  events  in  a  simulation  model. 

Many  VRTs  are  applicable  to  inventory  models.  Control 
variates,  antithetic  variates  and  conditional  Monte  Carlo 
methods  are  among  those  suitable  techniques  (6) . 

Ehrhardt  (6)  studied  the  use  of  each  of  those  three 
techniques  to  reduce  the  variance  in  estimating  the  para¬ 
meters  of  an  inventory  model.  He  combined  two  of  them  to 
achieve  more  reduction  in  variance  of  estimation.  He 
concluded  that  conditional  Monte  Carlo  is  the  best  sole 
technique  in  reducing  the  variance  when  applied  to  the 
underlying  inventory  model.  He  also  concluded  that  combin¬ 
ing  conditional  Monte  Carlo  with  either  of  the  two  other 
techniques  would  improve  the  variance  reduction  attained  by 
the  sole  technique.  Those  conclusions  were  drawn  from 
experimentation  with  the  following  inventory  model. 

The  Inventory  Model.  A  multi-item  inventory  system  which 
is  observable  at  discrete  intervals  of  time  was  studied. 

Each  of  its  items  has  an  did)  demand  from  one  period  to 
another.  An  order,  when  placed,  is  delivered  after  a  fixed 
number  of  time  periods  L  ,  and  any  unfulfilled  demand  in 
a  period  is  backlogged  to  be  satisfied  later. 
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The  inventory  cost  in  the  model  consists  of  a  fixed 
set  up  cost  per  order  k  ,  a  holding  cost  per  unit  per 
period  h  ,  and  a  penalty  cost  per  backlogged  unit  per 
period  P  . 

Only  partial  knowledge  of  the  demand  distribution  is 
available  to  the  decision  maker.  Two  policies,  containing 
different  degrees  of  information,  were  considered.  Each  of 
the  two  policies  is  of  (s,S)  type  where  s  is  the  lower 
inventory  replenishment  level,  and  S  is  the  upper  inventory 
level. 

The  first  policy,  called  the  Empirical  Normal  Approxi¬ 
mation  policy,  requires  the  knowledge  of  only  the  mean  and 
the  variance  of  the,  assumed,  normally  distributed  demand, 

M  and  o2  .  This  policy  can  also  be  called  the  constant 
policy  since  s  and  S  are  considered  constant  in  this 
case.  In  fact,  expected  values  of  the  operating  character¬ 
istics  of  the  system  can  be  calculated  directly  without 
simulation  in  this  case  using  an  analytical  approach,  but 
the  author  (6)  used  simulation  only  to  compare  this  policy 
with  the  second  policy. 

The  second  policy  considered  is  called  Statistical 
Normal  Approximation  policy.  In  this  case,  only  sample 
statistics  of  demand  are  available.  The  decision  maker  has 
to  revise  his  policy  periodically  since  he  would  not  know 
that  the  demand  distribution  is  stationary.  During  each 
revision  interval,  the  sample  mean  and  the  variance  of 
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demand  are  computed.  These  statistics  are  substituted  for 
the  actual  mean  and  variance  of  the  demand  distribution  to 
give  an  (s,S)  policy  to  be  used  in  the  next  interval  where 
new  statistics  are  collected.  It  was  required  to  develop 
efficient  simulation  techniques  to  evaluate  that  policy. 

The  simulation  experiment  was  designed  to  estimate  mean 
values  of  the  following  five  operating  characteristics  of  the 
system. 

1.  Holding  Quantity:  the  average  number  of  units 
on  hand  at  the  end  of  a  period. 

2.  Stockout  Quantity:  the  average  number  of  units 
backlogged  at  the  end  of  a  period. 

3.  Stockout  Frequency:  the  fraction  of  periods  in 
which  demand  is  backlogged. 

4.  Ordering  Frequency:  the  fraction  of  periods  in 
which  an  order  is  placed. 

5.  Total  Cost:  the  average  total  cost  per  period. 

The  objective  of  the  study  was  to  identify  the  suitable  VRTs 
that  yield  low  variance  estimate  of  the  expected  values  of 
the  five  operating  characteristics  for  a  given  cost  of 
computation. 

Simulation  Techniques.  To  estimate  the  expected  values 
of  the  five  operating  characteristics,  four  simulation 
techniques  were  used;  crude  Monte  Carlo  (direct  simulation), 
antithetic  sampling,  control  variates  and  conditional  Monte 
Carlo.  When  using  each  of  these  techniques,  the  vector 
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X  =  {x^  i=l,2,...,N}  represents  the  stochastic  process  of 
output  of  an  N-period  simulation  of  an  operating  character¬ 
istic.  For  any  of  the  five  operating  characteristics,  the 
expected  value  will  be 

1  N 

U  =  limit  k-  l  x.  (5.1) 

N-*-“  i=l 

Using  the  simulation  method,  j  will  estimate  jj  as 
for  N  periods. 

Direct  Simulation  (Method  1) .  Using  crude  Monte  Carlo 
the  actual  realizations  are  collected  as  simulated  and 
is  estimated  as 


1 

N 


N 


l 

i=l 


(5.2) 


the  variance  of  this  statistic  is 


(1)  ,  N  .  N  N 

Var  \i  -  —  l  Var  X.  +  —  l  l  .  Cov(X.,X.)  (5.3) 

N  N2  i=l  1  N2  i=l  j=l,J^l  1 


Antithetic  Sampling  (Method  2) .  To  apply  the  antithetic 
variates  in  this  study,  direct  simulation  was  applied  to  the 
first  N/2  periods,  then  the  simulation  was  restarted  where 
the  antithetic  variates  were  used  to  the  second  N/2  periods. 
That  is,  if  the  set  U  of  uniform  deviates,  CUi #  i*l, 2, . . . ,N/2} 
is  used  to  generate  the  first  N/2  demands,  then 
{ (1-U^) ,  i«l,2, . . . ,N/2)  is  the  set  of  deviates  used  to 
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generate  the  second  N/2  demands.  The  estimate  of 


is 


(2) 

V 

N 


I  X.  +  X 

i=l  1 


(5.4) 


where  is  the  operating  characteristic  output  when  using 

the  antithetic  stream  of  deviates.  The  variance  is  given  by 


(2)  2  N/2 

Var  y  =  —  £  [Var  x.  +  Cov(x.,x^)] 

N  N2  i=l  1  1  1 


2  N/2  N/2 

+  —  l  l  [Cov (x. ,x . )  +  Cov(x. ,x^)]  (5.5) 

N2  i*l  j=l  1  3  13 

i^j 

To  achieve  a  reduction  of  the  variance,  the  covariance  terms 
in  (5.5)  should  be  sufficiently  negative  to  make  Var'  1  < 

N 

Var^1^  Which  is  not  always  true. 

Control  Variates  (Methods  3a  and  3b) .  To  apply  control 
variates  approach,  an  approximate  model  should  be  simulated 
along  with  the  model  of  interest.  For  a  given  operating 
characteristic  x*  =  {xt,  i=l,...N)  denotes  the  stochastic 
process  of  output  for  the  approximate  model  simulation,  and 

l  N 

V*  -  Lin  i  l  x*  (5.6) 

N-*-«  w  i-1  1 

is  the  expected  value  of  this  output,  y*  is  assumed  to  be 
known  exactly.  The  control  variates  approach  estimates  y  as 
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(3) 

W 

N 


1  N 

£  l 

N  i=l 


(5.7) 


where 


ri  =  xi  "  6(xi  ~  (5.8) 


(2) 

and  B  is  a  constant.  The  variance  of  y  is  given  by 

N 

(3)  .  N 

Var  y  =  —  l  [Var  x.  +  B2  Var  x*  -  23  Cov(x.,x*)D 

N  N2  i=l  x  1  11 


+ 


,  N  N 

—I  l  [Cov(x.,x.)  +  B2  Cov(x*  x*) 

N2  i«l  j=l,j*L  1  3  1  3 


-  B  Cov(x|,Xj)  -  B  Covfx^x*)] 


(5.9) 


This  variance  depends  on  the  chosen  value  of  the  constant  p 

It  is  difficult  to  determine  the  value  of  3  that  will 
(3) 

minimize  Var  y  .  The  study  used  two  common  approaches 

N 

for  the  choice  of  B  •  Method  3a  uses  the  regression 
estimate  of  B*  while  B  is  set  to  be  one  in  Method  3b. 

In  both  cases,  the  expected  values  of  the  operating  charac¬ 
teristics  from  the  Empirical  Normal  Approximation  are  control 
variates  for  the  statistical  policy  simulation. 
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Conditional  Monte  Carlo  (Methods  4a,  4b,  4c  and  4d) . 
To  apply  the  conditional  Monte  Carlo  technique  to  the 
considered  inventory  model,  the  following  notations  have 
been  used.  Let  w  =  (w^,  i=l,...,N>  be  a  sequence  of 
vectors  specifying  the  state  of  the  simulation  at  each 
period.  Let  v  =  {v^,  i=l,...,N}  be  a  sequence  of  condi¬ 
tional  expectations 


v±  =  Efx^A^) 


(5.10) 


The  estimate  of  y  is  given  in  this  case  by 


(4) 


N 


N 


1  V 

»  Jr* 


(5.11) 


with  a  variance  of  estimation  given  by 


(4) 


Var  y 


N 


.  N  .  N  N 

l  Var  v.  +  -i-  l  l  Cov(v. ,v.) (5.12) 
N2  i=l  1  N2  i=l  1  3 


The  estimate  of  y  and  the  variance  given  by  (5.11)  and 
(5.12)  can  be  used  for  each  of  the  five  operating  character¬ 
istics  considering  the  following  notations.  Let  D  and 

be  the  demand  in  one  period  and  in  (L+l)  periods. 
Let  Y ^  denote  inventory  on  hand  and  on  order  in  period  i 
and  let  s^  be  the  value  of  reorder  policy  in  period  i  . 
Finally,  let  (a+)  be  the  max  of  a  and  0  .  Using  those 
notations,  the  conditional  expectation  for  the  holding 
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quantity  is 


V^1]  =  E[(Yi  -  D*(L+1))+D  (5.13) 

and  for  stockout  quantity 

vf2)  =  e[d*(L+1)  -  Y.)+3  (5.14) 

for  stockout  frequency 

V*3)  =  Pr  (d*(L+1)  >  Yi)  (5.15) 

for  ordering  frequency 

V^4)  ■  Pr  (D  >  (  Y±  -  s±))  (5.16) 


and  for  total  cost 


v<5) 

i 


hvf> 


+  P  v|2) 


+  k  v]4) 


(5.17) 


The  functions  given  by  (5.13)  through  (5.16)  were  calculated 
before  simulation  for  a  feasible  range  of  arguments  and  then 
used  appropriately  in  each  simulation  period  as  the  condi¬ 
tional  estimates. 

Four  variations  of  conditional  Monte  Carlo  were  examined 
in  the  study  (6) : 
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( 
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Method  4a: 


Method  4b: 


Method  4c: 


Method  4d: 


Conditional  statistics  are  collected  for  all 
operating  characteristics. 

Conditional  statistics  are  collected  for  all 
operating  characteristics  except  ordering 
frequency  which  is  simulated  directly. 
Conditional  statistics  are  collected  for 
operating  characteristics  except  holding 
quantity  which  is  simulated  directly. 
Conditional  statistics  are  collected  for 
only  stockout  quantity  and  stockout 
frequency. 


In  each  of  the  above  cases,  total  cost  is  taken  as  the 
weighted  sum  of  the  other  characteristics. 

Combined  Methods  (Methods  5a  and  5b) .  When  the  results 
of  Methods  1  through  4  were  available,  the  author  combined 
<_..tithetic  sampling  and  conditional  Monte  Carlo  to  get  two 
mixed  methods: 

Method  5a:  Method  4a  with  antithetic  sampling 

Method  5b:  Method  4b  with  antithetic  sampling 


Summary  of  the  Results 

When  simulation  was  performed  for  all  mentioned  methods, 
variance  reductions  gained  by  applying  different  VRTs  com¬ 
pared  with  crude  Monte  Carlo  technique  were  calculated. 
Antithetic  sampling  and  control  variates  were  found  to  yield 
meager  results,  while  conditional  Monte  Carlo  technique  gave 
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superior  performance.  It  reduced  the  variance  of  stockout 
quantity  estimates  by  factors  ranging  from  1.0  to  7.6  for 
the  statistical  policy  and  from  5.4  to  740  for  the  constant 
policy.  When  conditional  Monte  Carlo  was  combined  with 
antithetic  variates  method,  the  corresponding  variance  reduc¬ 
tion  ranges  are  1.4  to  12  and  3.8  to  460,  respectively.  The 
power  of  the  applied  VRTs  was  also  significant  for  estimates 
of  the  aggregate  total  cost.  Specifically,  conditional 
Monte  Carlo  reduced  the  variance  of  total  cost  by  factors 
ranging  from  2.1  to  14  for  the  statistical  policy  and  from 
2.3  to  480  for  constant  policy.  When  combined  with  antithetic 
sampling,  the  corresponding  variance  reduction  ranges  are 
4.8  to  20  and  2.0  to  1300,  respectively. 

The  author  pointed  out  that  the  cost  of  computation  was 
nearly  the  same  for  all  variance  reduction  schemes,  because 
the  computational  effort  was  dominated  by  updating  the  state 
of  the  system  and  by  output  analysis  in  this  study.  This 
means  that  the  above  given  variance  reduction  factors  can  be 
a  direct  measure  of  the  efficiency  of  corresponding  VRT 
relative  to  crude  Monte  Carlo. 

Application  of  VRTs  in  Queuing  Problems 

Another  example  of  the  fields  of  application  of  VRTs  is 
the  simulation  of  queuing  problems.  Many  studies  have  been 
accomplished  concerning  the  application  of  VRTs  in  this  field. 
One  of  these  studies  (10)  was  done  by  Gaver  where  different 
Monte  Carlo  techniques  were  discussed  and  then  applied  to 
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^  queuing  examples.  The  bases  for  straightforward  Monte  Carlo, 

antithetic  variables,  stratification,  control  variates,  and 
concomitant  variables  methods  were  summarized  in  the  paper 
The  basis  for  the  last  method  is  worth  mentioning  since  it 
has  not  been  discussed  before. 

Concomitant  Variables 

Suppose  that  realizations  of  the  random  variables 
X  (inter-arrival  times,  service  times,  etc.)  are  used  to 
create  the  realization  of  W  (waiting  time  of  an  individual 
in  the  queue) ,  where 

w(j)  .  f(x<j))  (5.18) 

Commonly,  and  are  monotonically  related  and 

then 

Cov  =  ci  (5.19) 

where  c.^  is  either  positive  or  negative.  In  fact, 

Edxf^D  =  ElIX,^  (5.20) 

since  X  is  a  given  specified  input.  When  sampling  only 
k  times,  the  realized  X  -values  will  deviate  from  their 
means.  Then  a  linear  correction  to  simple  average  will  be 
needed 

(  > 
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(5.21) 


1 

E[wDc  -  w  +  l  ya  (Xi  -  E[xiJ) 

where  can  be  estimated  in  terms  of  the  covariance  of 

W  and  X^  ,  and  the  resulting  estimate  is  unbiased  and 
consistent  asymptotically  as  the  sample  size  increases. 

Queuing  Examples 

The  VRTs  mentioned  before  were  well  illustrated  by 
consideration  of  a  very  simple  queuing  problem.  The  waiting 
time,  wr  ,  of  the  n-th  arrival  to  a  single-server  facility 
can  be  written  as 

w  =  max  [w„  t  -  A.  +  S.  , ,  0]  (5.22) 

n  n— x  n  n*“-L 

where  Ar  is  the  inter-arrival  period  elapsing  between  the 

(n-l)th  and  n-th  arrivals  to  the  queue,  and  Sn  is  the 

service  time  of  the  n-th  customer.  If  (A  }  and  {S  >  are 

n  n 

mutually  independent  sequences  of  (iid)  random  variables 

with  E[A  3  =  E[A]  >  E[S  3  =  E[S3  ,  and  if  other  moments 
n  n 

exist  as  required,  then  a  stationary  distribution  for  wr 
exists  as  n  -*•  «>  .  The  behavior  of  the  system  depends  on 

the  relation  between  ECAU  and  E[s3  ;  also,  it  depends  on 

the  number  of  arrivals  n  .  The  following  cases  were  dis¬ 
cussed  by  Gaver  (10) . 


€ 


V-13 


J 


1.  When  E[a3  is  little  larger  than  eCs3  ,  the 
queue  tends  to  be  long  developing  "heavy  traffic" 
situation.  In  this  case,  approximate  solution  of 
the  system  based  on  diffusion  equation  solution 
is  available. 

2.  When  ECA3  <  Ets3  the  queue  tends  to  grow  and 
little  information  is  available.  In  this  case, 
(5.22)  will  be  simply 


( 


w 


n 


=  w  .-A  +  S„  , 

n-1  n  n-1 


(5.23) 


where  w*  will  be  approximately  normal  if  Ar 

and  S  have  finite  variances.  The  mean  of  w„ 
n  n 

in  this  case  is 

ECw  1  =  (n-1) (ECS3  -  E[a3)  ,  w=  0  (5.24) 

n  j. 

For  small  and  moderate  values  of  n  ,  the  variances 

of  A  and  S  may  not  be  finite,  then  one  would 

have  to  use  simulation  to  estimate  ECw„3 

n 

The  author  (10)  applied  various  VRTs  to  estimate  ECwn3  for 
selected  values  of  n  ,  focusing  on  control  variates  and 
concomitant  variables  approaches.  The  following  numerical 
example  was  used  to  display  the  effect  of  selected  VRTs  on 
the  accuracy  of  estimation. 


;) 


$ 

4 
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Numerical  Example.  In  a  single  server  queuing  system. 


service  times  are  taken  to  be  exponentially  distributed 
(with  mean  =  10/q)  and  the  interarrival  times  of 

customers  are  taken  to  be  constant  (regular  arrivals)  or 
exponentially  distributed  (poisson  fashion  arrivals) ,  with 
unity  mean  in  both  cases. 

In  spite  of  the  apparent  simplicity  of  such  a  system, 
its  transient  response  is  not  easily  characterized  mathemat¬ 
ically.  Simulation  should  be  used  to  estimate  the  parameters 
of  the  system  or,  alternatively,  the  diffusion  approximation 
would  be  used.  Results  of  applying  the  selected  VRTs  and 
diffusion  approximation  are  tabulated  in  Tables  5.1  and  5.2 
for  the  regular  and  poisson  arrival  cases.  Discussion  of  the 
results  in  those  tables  is  given  below. 

Discussion  of  Results 

Rows  (1)  and  (2)  of  Tables  5.1  and  5.2  show  the 
results  obtained  when  25  independent  realizations  were 
averaged  to  estimate  E£wn3  .  The  same  random  numbers  were 
used  to  apply  antithetic  variates  approach  to  estimate 
wR(a)  and  its  variance,  rows  (3)  and  (4).  Comparison  of 
variances  in  rows  (2)  and  (4)  indicates  that  antithetic 
variates  approach  has  produced  an  improvement  even  after 
considering  the  additional  computational  effort  when  simu¬ 
lating  a  total  of  50  realizations.  The  improvement  in 
Table  5.2  is  smaller  than  that  in  Table  5.1  due  to  the  added 
variability  contributed  by  the  random  arrivals. 
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{  The  simple  control  variable  technique,  rows  (6)  and  (7), 

gave  better  improvement  for  large  values  of  n  than  anti¬ 
thetic  variates  approach.  This  appears  clearly  in  Table  5.2. 

Rows  (8)  and  (9)  display  the  effect  of  adjusting  the  straight¬ 
forward  estimate  in  accordance  with  the  concomitant  variable 
that  equals  the  sum  of  the  first  n  service  times  in 
Table  5.1,  while  Table  5.2  considers  both  service  and  arrival 
times  as  concomitant  variables.  The  behavior  of  the  concomitant 
variables  technique  was  similar  to  the  control  variables  and 
antithetic  variables.  Rows  (9)  and  (10)  display  the  results 
of  applying  concomitant  variables  to  the  components  of  the 
antithetic  estimate  of  (3)  and  (4) .  Results  in  this  case  are 
better  than  any  of  the  above  cases.  Rows  (12)  and  (13) 
indicate  the  value  of  regression-adjusted  control  procedure 
where  regression  was  used  to  determine  the  value  of  Bq 
used  in  the  modified  estimator  of  the  mean  waiting  time. 

E[wD  =  w  +  BQ  (w*  -  E[w*D)  (5.25) 

Row  (4)  shows  the  estimation  of  the  mean  of  w  when 

n 

diffusion  approximation  was  used.  The  results  in  row  (14) 
agree  quite  closely  with  control  variables,  row  (6),  and 
regression-adjusted  estimates,  rows  (8)  and  (10),  for  large 
n 
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Application  of  VRTs  in  Computer  System  Performance  Measurement 

Another  example  of  fields  of  applications  of  VRTs  is 
computer  system  performance  measurement.  A  study  in  this 
field  (7)  was  performed  where  the  antithetic  variates  approach 
was  employed  to  reduce  the  variance  of  estimation. 

Performance  evaluation  of  a  computer  system  is  very 
important  in  both  design  and  utilization  phases.  In  the 
second  case,  the  most  important  and  difficult  problem  in  the 
design  of  a  measurement  experiment  is  the  determination  of  its 
time  duration  since  the  increase  of  this  duration  increases 
the  cost  of  the  experiment.  On  the  other  hand,  the  time 
duration  of  the  measurement  experiment  should  be  long  enough 
to  yield  a  good  estimate  of  the  unknown  parameters.  This 
contradiction  is  similar  to  that  which  arises  in  the  simulation 
problems  where  the  sample  size  utilized  should  be  large  enough 
to  obtain  a  good  estimate  of  the  unknown  population  means 
and,  at  the  same  time,  this  sample  size  should  be  small  enough 
to  keep  the  simulation  cost  feasible. 

The  above  problem  was  solved  in  simulation  by  utilization 
of  the  suitable  VRTs  which,  when  applicable,  give  the  required 
precision  using  a  small  sample  size.  This  fact  encourages 
the  use  of  such  techniques  in  the  field  of  computer  performance 
measurement.  To  examine  the  profitability  of  applying  VRT 
in  this  field,  an  experiment  with  an  existing  synthetic  job 
generator  (described  later)  for  the  computer  system  was 
designed  and  performed  using  antithetic  variates  as  a  VRT. 
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The  Synthetic  Job  Generator 

The  synthetic  job  generator  used  in  this  experiment  was 
designed  to  model  the  I/O  (Input /Out put)  behavior  of  user 
jobs  in  a  CDC  6400  computer  system.  The  values  of  the 
parameters  of  this  generator  were  measured  by  employing  a 
trace  technique  which  records  a  complete  history  of  CPU 
(central  processing  unit)  containing  the  time  switched  between 
user  jobs  and  a  history  of  the  two  disks.  The  synthetic  job 
generator  was  designed  according  to  the  information  obtained 
from  tracing  the  re<tl  system  for  a  period  of  eight  hours  of 
normal  production  processing.  Then  the  designed  synthetic 
job  generator  was  used  to  perform  the  following  experiment. 

The  Experiment.  A  single  performance  measure  for  the 

system,  namely  the  mean  lob  elapsed  time  (i.e.,  the  mean 

time  a  job  generated  by  the  synthetic  job  generator  spends 

in  the  system)  was  used.  The  population  is  composed  of  the 

elapsed  times  of  all  possible  jobs  the  generator  may  produce. 

Its  mean  is  to  be  estimated  using  a  small  sample  size.  Let 

tR  be  the  mean  of  a  sample  of  size  n  .  If  many  samples 

of  the  same  size  were  drawn  from  the  population,  then  tR 

will  be  a  random  variable  with  a  mean  denoted  by  E(t  )  and 

n 

a  variance  denoted  by  Var(t  )  .  To  use  the  method  of 

n 

antithetic  variateB,  two  samples  of  N/2  jobs  each  are  used. 
Let  t^/2  be  the  mean  of  the  first,  and  t^j  be  the  mean 
of  the  second.  The  antithetic  estimator  of  the  population 
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-  °  tn/2  +  (1  ‘  “>  4n/2  •  0  <  o  <  1 

and  the  variance  of  t*  is  given  by 

Var (t*)  =  a2  Var(tn/2)  +  (1  -  a)2Var(t'/2) 

+  2a (1  =  a)  Cov(tn^2, t'^2) 


(5.26) 


(5.27) 


To  establish  the  negative  correlation  between  the  random 
variable  realizations  in  the  two  samples,  the  random  variable 
values  in  the  first  sample  is  generated  using  a  sequence  of 

random  numbers  ,  R2 . Rn/2  from  U(0,1)  and  those  of 

the  second  sample  are  generated  using  (1-R^) , (1-R2) , . . . (1-R^)  . 
In  this  experiment,  n  was  chosen  to  be  80  jobs,  40  jobs 
for  each  sample,  and  a  was  chosen  to  be  1/2.  When  the 
obtained  results  were  analyzed,  the  authors  drew  the  following 
conclusions. 

Conclusions.  The  authors  concluded  that  their  experiment 
to  demonstrate  the  feasibility  of  applying  antithetic  variates 
method  was  not  completely  successful,  but  they  also  concluded 
that  this  lack  of  success  was  due  to  the  instability  of  the 
devised  synthetic  job  generator  they  used.  It  is  the  opinion 
of  the  authors  that  stochastic  job  generators  can  be  designed 
which  have  all  the  properties  required  to  make  the  method 
of  antithetic  variables  very  effective  in  reducing  the 
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variance  of  an  estimator  and,  hence,  for  a  given  confidence 
interval  and  a  given  confidence  level,  in  decreasing  the 
minimum  duration  of  a  measurement  run.  The  authors  also 
stated  that  better  results  would  have  been  obtained  if  values 
less  than  1/2  had  been  chosen  for  a 

Comment .  It  is  clear  from  the  authors'  conclusions 
that  the  application  of  the  antithetic  variates  method  in 
their  experiment  was  not  completely  successful  for  the 
reasons  given  above.  It  would  be  better  if  they  tried  to 
use  other  values  of  a  to  minimize  the  variance  for  a  given 
status  of  the  experiment.  They  could  also  have  obtained 
better  results  if  their  synthetic  job  generator  was  more 
stable  and  had  the  features  suitable  for  applying  antithetic 
variates  approach.  Another  possibility  would  be  to  try  other 
applicable  VRTs. 


VI.  Conclusions  and  Recommendations 


VRTs  were  investigated,  clarified,  contrasted  and 
illustrated  in  the  former  chapters.  From  the  extensive 
literature  review  and  the  computed  numerical  examples,  it 
was  found  that  none  of  the  known  VRTs  is  generally  superior 
to  the  others.  In  other  words,  there  is  no  single  technique 
which  is  the  most  suitable  technique  for  every  simulation 
problem.  The  condensed  table  (Table  3.2)  of  VRTs'  charac¬ 
teristics  identified  the  fields  of  application  of  each 
technique.  It  is  clear  that  various  techniques  can  be 
applied  to  the  same  problem  in  certain  fields.  In  this  case, 
one  can  compare  the  efficiency  of  various  techniques.  The 
most  efficient  technique  is  the  one  which  utilizes  most 
of  the  available  information  about  the  underlying  simulation 
process  and  gives  the  most  accurate  estimation  of  the  para¬ 
meters  (minimizes  the  sample  variance)  with  minimum  computa¬ 
tional  effort.  In  many  cases,  these  three  extreme  objectives 
cannot  be  achieved  simultaneously  by  one  of  the  VRTs.  Only 
the  involved  analyst  can  weigh  them  to  choose  the  optimal 
technique  suitable  for  his  simulation  problem,  objectives, 
and  available  resources.  It  is  not  practical  that  an  analyst, 
handling  a  real  world  problem,  will  apply  several  applicable 
VRTs  to  his  problem  and  then  compare  them  to  choose  the 
optimal  technique.  A  study  like  this  and  other  studies 
devoted  to  the  application  of  VRTs  to  the  specific  kinds  of 


problems  can  help  the  analyst  select  a  profitable  VRT.  The 
degree  of  profitability  will  then  depend  on  how  successfully 
the  technique  is  implemented. 

A  possible  useful  extension  of  this  work  is  an  organized 
collection  of  real  world  problems  in  different  fields  of 
applications. 
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Computer  Program  Listing  with  the 
Results  for  Russian  Roulette  Example,  Case  1 
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NO.  OP  REDUCTION?  --  366  F ( RS  GT  PX  )  -  .25CC 

*  v  *  4  A  *  *  *•  *  4  *  4  4  4  4  4  *'  4  4  *  *  *.  *  4. 4  4. 4:  *  4  4  *  *  4  4  4  *  r  4  4  4  *  *•  4  *  *'  4  4  >.  4  4 

PX  =  .550  FT  r  .700  FT  -  .800  P3“ . 800 

5  *  ;  *  4  4  .  *  *  *  *■  *  *•  *■  v  > ;.  /  *  x  *■  *  *  * *  *  *  *  *  4  *.  *  4  „  *  *  *  *  *■  *  *  >•  *  *  >;  ■  *.  *  « . 

NO.  Or  RE  DUCT  J  ONE  -•  5  <;2  PiO  (  T  r  X  )  .1-55: 

*4**.**£*.*  ».****444  44*  r***x*i  a*;:4;;  ;  '  .*•*:*  *•*  **•*  s  ****•»> 

PX  -  .  60' ■  FT  *  .700  FT  -  .800  P3  =  .  800 

*  *•*  *******  *******  *  i  44*  *  *.*  *4*4  4**4**  4  ■  i  4  44  :>  *  4  *»  *»  4  > 


B-3 


F‘  ( R'S  C-T  py  > 


er.30 


NO.  C*r  REI-.ICT3  ONS  =  623 

*4>  -  4**4*  *4*:?.  *4*  44 4*  *>  *.*  ***** * ****  *  i  **»•* 

F'X  -  ,650  F*:  -■  •?_••'■  F*2  ~  ,  Pi'-.'  r  vr.  .300 

**  4*4  O'*  4*44*  *44*  4  y  :.  -•  *  *4  ■'. v '  *  *•#  *  ».v:>  *.*  *4*4,  » ,j  4./  ■•■  * 

NO.  OF  i'  EDL'C!  IONL  -  %122  F<RS  G"i  P‘>  )  =  ,0260 

4  4. 4  4  **X  *  i  S  4  4  4  4  4  *  4  *.t  4  >.4'  *  *4444*4  *  4  4  X'X  *  X  y  X  ,f  *  V  '  A'l)<  4  >-.  Jl 

<  ”  4<  4'  X  4  *4  4  *  4  4  *44  4  4  *  *  Js  *  *  4  4  4. 4  *  *  >  4  4444  *44  *  4  *  4  v-  v  A '  4  >■  *  * 

no.  (.■“  :: oi o: -ion?  -  r--,3v  r<n<:  or  py  '  =  , oooo 

**<  4  44  4444  4,f  *  4  44  x  y  *  4.444  t  44  >  4  4  4  44*4  4  ::.**  4  v  .4  4  »  4;  t  4  i  4 

py  -  r  •.  ■  .  ~>.n*  f  :  -  .300  r-  ;•••  . oc 

4  4  *  *.  *  X  4  4  *  4  4  *  j»  *  y  4  *  4  »  4  4  *  4  *  4  4  *  *  *  •*  »  *  *:  X  4  /  x  4  *'  4  1  ■  ,■.  X  X  4  •>  I 

NO.  OF  FiF.Iil.lC~  It*N6  *•  it  54  f  •  i.5  5T  P>.  •■•  ,  C 0  20 

4: 4  *4*4 4. 4* x.4 4 *4**4  4  4  4.444444  •.  t  *  *  4>  *  4  *  *’'  ‘  Ox  4  *  /  *  •  .;•  y  *  x 

i  r.n.zii 

'<>(»  ?■' 0x1  '-'Jo  Evt-ctr  i •::■■;  . 

'22  c.r  zr  cc  ri:  txFCiHjn  '•  . 


APPENDIX  C 


Computer  Program  Listing  with  the 
Results  for  Russian  Roulette  Example,  Case  2 


100“ 

'*r:c:.:.F- 1  7: 

:  io=-e 

i  » yy,  *  y  y  y  >•  >  >  ,*  no  i 

>*.**>  ■  (  1  1  i  -  y  ;■  v  , 

.  i  -  -  J  .  .  ■;  •  v . :  . 

-  *  ."1,  u  u.  '  -  f  , 

:i  30~  c 

Thu  rvoTtir  r  j  2: r  ~  o 

OF  7  7  i  F  ;!  LD  jii-1 

1  4  (>  ;;; 

>{*  f.  V.  Jf  %  *  >;  .  ill  •  V  .  •<  ;  v  y 

X-  L  •"  a.  H  /•  lit  j  *  .  v  , 

r.'l'Y  L  ■■ i  ' 

5: 7  0“ 

.1.  ■. 

•  V:4  •'r'.  C-  '■  ,  I  •  :  ,2. 

:i  •  •  o  •" 

r  x>* ,  <v 

:  ;  ()  f  \ 

wi;;i  •;  y-  ,  i  5  ;• 

2:i0*' 

PMKT  2  ,  Ti-iHE  r  50 1 

BVBTL'r'r  ]N  f'CuAL  L  F  ' 

v  - 

LD  200  L . - i  ,  10 

330™ 

PX=F*X+  .  DOS 

">4 

T.-O 

. 

*1-  p  •  •  ..  •  *v  .•% 

*-  '  - 

260- 

WRITE  <*,  25  )F’V  ,  'T  7  (  7 

)  ,  J  •  :  v  3  ) 

n  \ ... 

A.  ..  •  ’ 

25  rCrifirY.  <  / , 10>. »  i ’■  X  ~ 

jf  1“  Av  «■  w  •  :r-‘  *  & 

“3=  ‘ 

2  f  >  o  sr 

4  .34.3) 

290“ 

DO  10  T - 1 , 3 

300- 

1 0  on. )-pi  ( i)*(  1  -r- 1.  •;  1  > 

'  /  r } 

.  ’•?  *•  o :::: 

DO  100  1*1,1000 

320™ 

CAL  L..  NORM  (  P,N ) 

J  3  0 

Xl=  r"l(lHI;i«SQRT(V( 

1  > ) 

340- 

1 F"  ( X 1 » GT  ♦  PX  )  T HE N 

3  SO- 

I  COUNT- 1  COUNT  ■{  1 

V.’  (.)  v' 

NF:  DC*  HR  DC  4  2 

370- 

GO  TO  100 

3S0* 

ELBE 

390- 

COLL  NORM ( PN> 

400“ 

X  2  ~  Pi  (2)  +  R  k  *  5  0  F;  T  <  M  < 

2 )  ? 

410“ 

Rl2*l-  <  l-Xl  >>•:(  :i-X2) 

420- 

I F ( R 1 2 ♦ 6 T « P X )  THEN 

430- 

1  COUNT* ICOUNT  +  1 

440- 

NRDC-NRDC ; 1 

450- 

CO  TO  100 

440- 

«...  i...  1  r 

4?'C,r- 

C ALL  NORM ( RN > 

‘■'0- 

X5-P1 ( 3 ) +HNV SORT < 0 ( 

3  >  > 

£  Q  0  • 

i  r  -  i-  ( :■  -xi )  x  >:  i-->:r>» 

.1  ■■■  X3  > 

soo  - 

IE  <  R5  .  GT  ,  PX )  '"HEN 

5 10- 

3  C0!|!\"T  “  irOiJKT  •!  1 

{  ;  j  *  \  ! 

1.  ?  <  *  •  1 

:  ,  ‘V,  0 

i  nd  ::t 

540- 

END  J : 

f  ,  ■ 

00  co:  ’tipt 

540“ 

Cl  *  j  COUNT.  -M  000.0 

C-l 


A  •  4  > 

61  0- 


WF.iTr-:*,j5> 

1 1>  CDftKAT  v  ■•,10X,S0( '*  •  • 

WRITE  Ot  ,20)  NRDC  *  C I 

20  r0'*M<VK3C*Xt  'NO.  O'  ftEDUL 710^8  ■■■  '  .IE  ,  3  v ,  ' 

-00  WRITE,  ,15) 

END 

SUI- ROUT T NE  N0Rh<F:N) 

A  =  0.0 

Cl=.02r'C99  77  A 

r*  ■  ■  /*. ,*\ ots.c; c.  /  c 

C3* • 0765429 1 2 
C 4 -  .257-408784 
r  *,  =  3,94  o  9  4  o  130 
DO  4  1*1,12 
CALL  RAND < ft; 

4  A  =  i*.  -f  ft 

ft; X  =;  (  A-6  .  )  /4 
R2~ftX*RX 

RN= ( < « <  C 1 *R2+C2  >  *R2+C3  >  *R2+C4 ) KR2+C5 ) *RX 

RETURN 

END 

SUBROUTINE  RAND (ft) 
ft CAL ft , S , XM 
INTEGER  IS 

DATA  S , XH/ . 23978 , 823 ♦ 53478/ 

S=S*XM 
IS»S 
R  =■':>- IS 

S'-M  .8234817 
RETURN 
END 


iN.rrrs 

63300  CM  STORAGE  USED. 

0.175’  CP  SEC  ONUS  COMPILATION  TIME, 

**XX##****»*##X*XXXX*XX*X»  *  4  X  4  X  *  X  *  X*  X'X  4  X  4  *  >*  XXX  X«  X* 
THREE  SJBBY STEMS  IN  PARALLEL 

PX  ~  .905  PI  "  .700  P?  ~  .700  P7  = .700 

X 4  4  # 5fC 4  X X X.  X  * X  *'  XX*X  X  X  X  X  A  if' A  A  X. j  X  A  >  A  %  X -1  >  <t  *  X  4  4-  *  *  A  X  if.  X  if-  A .»  4' 

NO.  OF  REDUCTIONS  »  571  F(F:E  GT  PX  i  •=  1 .0000 

X  4-4-4  4  t  ****  *  *  A  XXX  *  7  X  X  A  A  X  X  A  V  A  *  7  A  A  X .  7  X  *  4  X  *  4  7  *  4  X:  X  A  A  *.  *  X'  X 

PX  =  .910  PI  ~  .700  P2  -  .700  P3-.700 

4  *  X  *  XX*  A  X  A  *  A  4  *  4  7  7  4-  Si  4- 4  4-  %  4  4  X  4  4  X  X  A  X  *  A  A  X  X :  V  X  X  X  X  *  X"  *  *  X 4  4  4 
NO.  CF  PEI  uni  DNS  =  555  P(ES  GT  PX  )  =  1.0000 

4 4 4  4  4 4 4 X  *.  4  4 7  4 4 4 X: 444 4  4 4  4  4 X  4 4 X  X 4 4 4 4  4 4 4 4  X 4 4  4 4. X  4  4 X 4:  i  *  i 

PX  -  .915  PI  =•  <700  P2  *  .700  P3*.70C 

X 4 *  4 4 4 4  4  4 4  4 444-4 4- 4- X X X X X X X X X  X X X X X  X X X X X  X  X  X X X X  X XX*  X X  X. X 
NO.  OF  REDUCTIONS  *  497  P(RS  GT  PX  )  *  ,9990 

X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  4  X  *  X  X  X  X  X  X  4  X  X  X  X  X  XXX  X  X  X  4  X  X  >.  X  X  X 

PX  -  .920  PI  a  .700  P2  a  .700  P3« . 700 

*  X  4  4  X  X  X  X  4  X  X  X  X  X  x  4  X  X  X  4  X  X  4  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  x 
NO.  OF  REDUCTIONS  a  416  FCRS  GT  PX  )  =  .9980 

X X  X  4 X  X  X  X  X  X  X  X  X X  X  X  X  X X XXX  X 4 4 4 X  X  X X  X X  X X  X  X  X  X  XXX  X XXX  X  X  X X X 

PX  a  ,925  PI  •--  ,700  P2  a  ,700  P3«.700 

X  X X  X  X  X 4 X  X X  X  X  X X  X  X  X  X  X  X  X  X  X X  X X  X  X  X X X X  X X  X  XXX* X  X  X  X  X X  X  X X  X  4 
NO.  or  REDUCTIONS  a  394  P(RS  GT  PX  )  =  .9950 

X  X  X  X  X  x X  X  X  X  X  X  X  X  X  X  X  X  *  4  X  XX  X  X  XXX'.  XX  X  X  X  X.  XX  X  X  X  X X  X X  X  X  *  X X  X. x 

PV  ~  .  970  Pi  a  .700  P2  e  .700  P3«.700 

X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  *  X  X  X  4  4  4  X  X  X  X  X  X  X  X  X  X  X 

N'T.  OF  ::TD  ITT  JONG  ~  298  P(RG  GT  PX  )  =  ,9920 

XXX  X  XX  X  XXXXXXXXXX  XXX  X  XXX XXX  X  XX XX  v  XL X  "•*:  !  ;•  *:,;.*  *  '"X  S  4  44 
PX  ,935  PI  a  ,700  r?  =•:  .700  P3~  ,  700 


C-3 


J 


XXXXXX %  jfcjc jfi  MW#**  t  XXXXXX'X'XXXXXXXXXXXXX*  *  *  *.  X  X  X  *  v  *  **  *• 
NO.  OF  RET  DUCT  1  DNS  ~  25?  P(RS  GT  PX  >  -  .9820 

X  *  #  X  *  X  *  A  X*#*  X  XxXXX  X  *  X  X  *  X  X * X  x*#  X  *  *  *XX*  *  *  XXXX  x  sf  *  X  if.  A  X 

PX  -  .940  Pi  -  .700  P2  =  .700  P?=.700 

X  a:  X  X  X  X  *•  *  X  X  X  *  X  X-  X' v  •  X  '  •  X  XX**:*.  X  X  *  X'  X  X  *  X  X  *  X  X  X  -X  X  X'  *  X  *  X  x  X  X  X  x 
NO.  OF  REDUCTIONS  =  209  P(RS  GT  PX  >  -'  .9770 

XXXX*»XX*  X.  X  X  X  »:  X  *  X:  X  *  *  X  :r  X  X  X  *  X  X  X  X  #  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  V 

PX  =  .945  PI  «  .700  P2  ~  .700  P3=.700 

X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  x  X  X  X  X  X  X  X  x  X  X  X  X  X  X.  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X 
wo.  OF  REDUCTIONS  *  151  P(RS  GT  PX  )  ~  .9710 

X  X  X  x  X  X  X  X  X  X  X  X  x  X  X  X  X  X  X  x  x  X  X  X  X  X  X'  X  X  X  X  x.  X  x  X  X  X  X  X  X  X  X  X  X  X  X  X  *  X  X 

PX  »  .950  PI  *  .700  P2  =  .700  P3= . 700 

X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  x  X  X  X  X  X  X  X  X  X  X  X 
NO.  OF  REDUCTIONS  «  106  P(R<?  GT  PX  >  =  .9360 

X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X X  X  X  X  X  X  X X X  X  X  X  X  X  X X  X  X X X  X  X  X  X  X  XXX  X 
EK’D  RELID 

1.6000  MAXIMUM  EXECUTION  FL.. . 

4,323  CP  SECONDS  EXECUTION  TIME. 


C-4 


1.. ,  ft 


APPENDIX  D 

Computer  Program  Listing  with  the 
Results  for  the  Examples  of 
Systematic  Sampling  and  Stratified  Sampling 


<  oo~ 

1C  i=T 
1  10--C 

i20=c 

x  30-  C 
1 4  0  =  C 


1 70=  C 

ieo=c: 

1 90-C 

1 9  5  =c 
200  = 

:  10* 
220  = 

240= 

250= 

2  <50= 
270= 
200= 
290= 

20  0  = 
310= 
320= 
330= 
340  = 
350= 
360  = 
370  = 
380  = 
390= 
400  = 
410  = 
420  = 
430  = 


PROGRAM  MASTER 

*  *  *  *  *  *******  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  **********  *  *  *  *  ****  *  *  *  *  *  *  *  *  *  * 
THIS  PROGRAM  IS  AN  EXAMPLE  GF  USING, 

S Y  S T F: MAT  I C  S A rl P L I N 3  , 

STRATIFIED  SAMPLING  AS  V.r.T'S. 

IT  ALSO  PERFORMS  LRUIE  MG NTS  CARL  0  TECHNIQUE  TO  COMPARE 
THE  RESULTS  OF  THE  T HREF  METHODS. 

*  *  *  *  *  *  *  *  *  *  *  »•  *  *  *  *  *  ****  *  %  *  k  *  x.  *****  r  *  *  *  ******  *  *  *  *  *  *  *  *  *  *  * 

NOTE  : 

METHOD  =  1  CRUDE  MONTE  CARLO 

2  SYSTEMATIC  SAMPLING 

3  STRA T : F I ED  S A K P L I N G 
I NTEGER  I , J , K , N , NS <10), ME THOD 

R  E  ft  L  A  s.  <  1 0 > ,  A  h  <  1 0  > ,  ft  C  <  1 0  >,X,  Y ,  F  7 » 7  E  R  M ,  S  M  F ,  D  S  ,  F  U  N  C 

DATA  N  >  N , METHOD  756,4,2/ 

F  M = 0 , 0 
DO  1  1=1,50 
1  CALL  RAND  (R; 
ir<K,EQ.l >  THEN 
WRITE<* ,3) 

3  FORMAT <20X,  'CRUDE  MONTE  CARl.O  ■'  > 

NS  ( 1 )  =-N 
GO  TO  10 
END  IF 

DATA  <NS< J) , J=1 ,4)  74,4,4,4/ 

IF < METHOD, EQ. 2) THEN 
L'RITE  <  * ,  5 ) 

5  FORMAT < 20X, '  SYSTEMATIC  SAMPLING  ') 

GO  TO  10 

ELSE 

L’RITE<*,6) 

6  FORMAT < 20X , 'STRATIFIED  SAMPLING') 

DAT A < AD < J) , J= 1,4)  /0»0, ,36, ,62, . 83/ 

DATA < AL < J ) , J= 1,4)  /, 36 , ,62, ,83,1, 0/ 

GO  TO  110 

END  IF 


1 

440= 

10 

DS=1 ,0/K 

430  = 

DO  100  J  = 1 , K 

4  60= 

AL < J)=J*DS 

470  = 

100 

AD  <  J  >  =AL.  <  2 )  -  DS 

480= 

no 

DO  120  J = 1 , K 

A 

1.20 

ACi J)=AL( J)-AD<J) 

1  : 

500  = 

L’F  ITEOt  ,  130; 

1! 

510= 

130 

FORMAT ( 3X , 60 ('*')) 

)  i 

!  i  / '« 

520= 

UF-.I  TE  <  *  ,140) 

!i  ° 

530  = 

F  M ' ) 

140 

FORMAT C6X, ' J ' ,3X, 'I 

'  , 


D-l 


540- 

Ur;jTEOr,130> 

55  0  = 

SKP=0.0 

560” 

SH 1=0.0 

570« 

SM2=0 , 0 

560==' 

DO  205  J==3  ,  K 

590” 

SK..1  =0 , 0 

600 ~ 

HO  200  1=1  , NS  >'  J ) 

610- 

150 

CALL  RAND  C T  r.' 

620^ 

I F  ( RN .  GT  ,  tV_.  <  2 )  ,  OR .  RN .  Lt. .  AD  (  J )  .'= 

630  a 

3(3  TO  3  50 

640  = 

END  IF 

650== 

X*RN 

660= 

Y«AB(JHAC<J>*X 

670= 

FY=FUNC C  Y : 

660a 

SMi«sm+rr 

690  = 

FY2=FY**2 

700  = 

UR 1 7 E ■: * ,  1 6 0 > 2 , 1 , X ,Y,f Y 

710  = 

160 

FORMAT  ( 3 Y.  , 2 1 4 , 3 ( r  8 . 6 , 3X )  > 

y  o  /> 

200 

SM2*f.M2+FY2 

730a 

TERM*=AC ( J ) *8M1/N& ( J ) 

740=' 

FM=FM+TERM 

750a 

L'RITE(*,201)TERM,FM 

760== 

201 

r  OR  MAT  ( 43X  ,  2  <  F8 . 6 , 3X.  > ) 

770= 

205 

TERM 2=  < SM2*AC  <  J ) ( J)*-  TERM** 

700= 

F‘AR=FM 

790= 

VAR=TERM2 

000= 

WRITE ( * » 130) 

010  = 

WR I TL  < * ,  220 )  PAR , VAR 

620= 

220 

FORMAT <3X, ' EOT I MATE  OF  PAR  = ',28 

030= 

L'RITEU,  130) 

840= 

END 

850= 

FUNCTION  FUNC(Y) 

860  = 

REAL  FY* Y 

070= 

rUNO==  ( EXP  (  V  )-l. 03/1.718 

880= 

FY«FUNC 

890  = 

RETURN 

900= 

END 

910= 

SUBROUTINE  RAND(R) 

920= 

REAL  R,S,XM 

930= 

INTEGER  IS 

940== 

DATA  S r XM  7,23978, 823. 53478/ 

950= 

S=S*XM 

960= 

16=5 

970= 

R--B-1S 

980  = 

S=R» .8234617 

V90= 

RETURN 

1000  = 

END 

<  J  > ) 


'  varianc: 


D-2 


♦  * 


Jl 


/ 1 ;  |  ,  3  /  / 1  it  !  y  j  /  *  A 
1  r :  iAN  Cl.  5  ) 

run  .  file 

6  1  3  o  0  c  H  f  T  0  R  A  G  E  U  S  E  D . 

0.194  Cr'  SECONDS  COMPILATION  TIME. 


CRUEL  MONTE  CARLO 

51  *  *  *  x,  XXXXX*  *  X  **x  *  **»*#*  *  x:  »’  X  *  X  X  X  X  X  X  X  *  X  X  *  *  *  *■  *  *  #•  X  *  *  *  *  X:  *  X  X  X-  X  *  X  *  * 
.!  ;j  v  Y  FY  TERA  FM 

X X  X  X  *■  X * X  X X X: X  X  x  X  *  X X-  X  X X X: X  X  X X  * X  * X  X: X: X X X X X X X X X X X X  X X X X X X  XXX  X  XXX  X X * 


1 

.1. 

.20°<i CO 

. 209920 

.  135960 

j 

.  01?  5  4 13 

,  025423 

.014982 

3 

,  085874 

. 085874 

.052x9* 

.1 

4 

.8*9429 

.  86947° 

.  806489 

■j 

.  1  54523. 

,  1  5451  1 

.  097264 

■j 

6 

.  40292(1 

.  402926 

naacnn 

♦  «.  O  «...  *. 

j 

-■ 

.  97263* 

.  97263* 

.  957451 

8 

. 147175 

. 147175 

.  092292 

1 

. 353352 

. 353352 

.246701 

.i 

10 

, 146870 

,  .1  46870 

. 092086 

j 

•i  j 

.301687 

. 10167“ 

. 062309 

.1 

1 2 

,9  oo  or? 

. 900077 

.849632 

J 

■■ 

. 352683 

ro 

CD 

U  ' 

.246146 

1 

14 

. 595825 

.595825 

. 474113 

1 

15 

.831571 

♦  D  i  J.  /  4. 

. 754904 

h 

3  6 

.977236 

. 977236 

.964552 

.363494  .383494 

XXXXXXX* XXXXX  XX* XXXXXXX;#X XXX*  XXXXXX XXX XXXXX XXX* XXX* XXXX  xxxxx 

ESTIMATE  0“  FAR  .383494  VAR  3  ANCF  »  .119850 

X  X* XXXX  X  X  X  X  X X X X X X X  X XX X X X xxx X X X X XX X X X XXX*  *  xxxxx* y  XXXXX*  xXXX'X  y 
END  MABTEF: 

13300  MAXIMUM  EXECUTION  FL . 

0.043  CP  SECONDS  EXECUTION  TIME. 

ILF  QUOTA  EXCEEDED 


D-3 


RUN » FTN5 

6  ]”0'-  C“!  S7GRAGL  U=  El- . 

0 , 1  ? 5  C.  *  SEC 0 N It S  C D M f '  1 L A T 1 0 1 N  T I * "  , 

SYSTEMS"!  ]  f  SOM  PL  INS 

*4  *  **».#  »  *  **>****'#  *  *:*>  ** ******  *  **********  *  4 

j  i  >:  y  ry  term  fm 

%**)»  ****  ********  ******  **  *  *******  *****************  **<  ******** 


1 

j. 

.209920 

.052480 

.031363 

.025423 

.006356 

.003711 

1 

3 

.085879 

.  021 4 <;> B 

,012631 

4 

,154521 

.032630 

f\  7  r-  £ 

, 0044  1 4 

. 00 4 Hi  4 

*:> 

1 

. 402926 

,350731 

*■>  a  *  ir  "7 
•  *•-  *v  •'"»  v.‘  %;>  *. 

“> 

, 353352 

.  338338 

.  234351 

I' 

7 

♦  -  w  6  8  9 

.338171 

.234214 

r- 

4 

.38668.1 

. 346670 

.  24 1 1  32 

.059642 

.064057 

j 

.  594,19  3 

.  6  4  8  6 1 0 

.531362 

3 

.691941 

,  672986 

. 553836 

*7. 

7 

.537404 

.634351 

.515598 

*3 

4 

.720207 

,  680072 

, 566950 

, 135797 

.  199853 

4 

1 

.  938755 

,93468? 

.976122 

4 

2 

.806858 

.951715 

, 9 255 SO 

4 

"2 

,901946 

. 975486 

.961849 

j 

4 

.933105 

.983276 

. 973923 

.  239842 

,  439696 

*  *  *  *  4'  4  *  *  *  *■  4  *  *  *  V * 4  i  ****  *  4  *  *  *  4  4  4  *  4  4  4  '*:  4  4  4  4  4  4  4  4  4  4  4  4  >  4  4  4  4  4  4  4  4  4  4  4  4 
ESTIMATE  Dr-:  PAR  =  .439696  VARIANCE  -  .02207G 

44444  44444  44*  4  44444444  4  4 44*44444444 4  4  4  4 4 4  4  4  4  4  4  4  4  4  4  4  4  4  4  4 4  4  4  4 4 
END  MASTER 

13300  MAXIMUM  EXECUTION  FL. 

0.050  CP  SECONDS  EXECUTION  TIME. 
riLi:  QUOTA  EXCEEDED 


D-4 


*r 


RUN  r  FTN5 

61300  CM  STORAGE  USED, 

0.178  CP  SECONDS  COMPILATION  TIME. 

SIR  A  T I  r  i  r  r.  SAMP  L  TNG 

**#Y  **  YYYYYYYYY  ******'##•*  *********************  Y*Y  #*#**#.*****.*■* 

J  I  X  Y  FY  TERM  Fr; 


J 

I 

x 

V 

FY 

Y  Y  *  Y 

: 

*:  Y  Y  Y  Y  Y  Y  Y  *  Y  Y  Y  Y  Y  Y  Y  Y  Y  Y  Y 

Y YYYY YYYY Y S 

1 

•» 

j. 

. 209970 

.075571 

. 045693 

1 

.025423 

. 009152 

.005352 

3. 

-jr 

OS 5874 

.  03091  5 

.018276 

1 

4 

.  5  54521 

. 055628 

.033297 

1 

.402926 

.464761 

.344373 

'■j 

♦595825 

.514914 

.392023 

*v 

3 

.386681 

.  4  6 0 5 3  7 

.340468 

*•! 

4 

.594443 

.514555 

♦  J:  "/  1.  Cj  /  *1 

3 

i 

.691942 

.  669187 

3 

p 

. 755606 

. 778677 

, 686024 

w 

3 

.  822945 

.  792819 

, 704084 

3 

4 

.720287 

.771260 

. 676653 

4 

1 

.938755 

. 989588 

.983 7 7 6 

4 

n 

.901946 

.  98333.1 

. 974008 

4 

3 

.933105 

. 988628 

. 982273 

4 

4 

.962816 

. 993679 

.990194 

. 009235 


.0  95455 


.1  43637 


vn'y2/*r< 


1 0 4  6 90 


♦248327 


.167036  .415363 

Y:  6'  Y  Y  *  *:  #YY’  *  #  *  *  #  Y  Y  Y  *'■  Y  Y  Y  *  Y:  f  *  *  *  *  Y  Y  Y  *  #  f  *  Y  *  Y  Y  Y  Y  Y  Y  *  *  Y * Y  Y  Y  *  Y  f  Y  f  Y  *  Y  v  Y. 
ESTIMATE  OF  PAR  --  .415363  VARIANCE  =  .017465 

YYYYYY  *##*  YY  Y'YYYYY.YYYYYYYY  YYYYYYY  YYYY>  Y  Y  Y  s*  YYYY  :.YY  YY YYYYYY*  Y.Y 
END  MASTER 

13300  MAXIMUM  EXECUTION  FL. 

0.044  CP  SECONDS  EXECUTION  TIME. 

FILE  QUOTA  EXCEEDED 


APPENDIX  E 


Computer  Program  Listing  with  the 
Results  for  the  Example  of  Control  Variates 


» A 


j  00= 

j.  *•!.  v'  "■  L 

130=0 
140  = 

1  50= 

.1  60” 
170" 

1 80= 

.1  90« 

::oo~ 

230= 

2  4  0 = 

250“ 
260= 
270= 
280= 
290  = 
300= 
31  0= 
320  = 
330  = 
340= 
350= 
360= 
370= 
380  = 
390  = 
400= 
410= 
420= 
430  = 
440= 
450  = 
460= 

4  70  = 
480  = 
490= 
500  = 
510= 
520  = 
530- 
540  = 
550  = 
560= 


10  F 


50 


100 


1  09 


1. 1 0 


PROG 5 OH  C0KTF;’. 

X  i  X  XX  XX XX  £ X'  *’jj  x  X  £:£  X  X  XXX.XXXX  X  X X  XX  XX  XX  XX  X  X  X  X.  X  >  X  *  X X  X  X  X  X  X  X  XXX*  X 

t ! ;  I 3  P  R  0  5  H  A  . :  U  5  e:  s  c  o  n  t  r  c-  L  V  A  R  .1  A  T  e:  6. 5  V .  h  .  T  . 

XX  XX  XXXXX  XX  XXX  x  xxxxxx  XX  X  X X; X X  X X X  X X  X  X  X  X  xxx  X  XXX  X  *  X  X  *  xxxxxx  X  X 

INTEGE  R  1  , 

R  E  A  L  T  E  F:  5 .1  ,  T  E  R  M  2 ,  £  NT  1,8  H  7  2 ,  F  A  X ,  F  UN  C 
V  A  T  A  N  v  F 1  /  2  0  *  .  5  / 

SMT 1=0.0 
SMT2=0«0 
SMT3=0  *  0 
WRITE  (X<,  10) 

OF: HAT  (  7X » 5  0  (  •'  * '  >  > 
wf:i  te:  ■:  x ,  20  > 

FORMAT <  lt-X  *  ‘ R'\ '  ,i2X,  /  TERM  'TERM2  ,3X,  '  TERM 3  '  I 

WRITE  OK, 10) 

ro  100  i* i,n 

CALL  RANIi(RN) 

TERM 1* FUND <RN> 

SMTlaSKT :l  4 TERM1 
TERK2=FAX<RN> 

TERM3= < TERM! -TERM2 ) XX2 

WRITE  <*,50) RN , TERM) , TERM2 , TERMS 

FORMAT (5X» 4 <4X »F9 . 6 > ) 

SMT3-SMT24 TERM3 
£  M  T  2 = S M  T  2  4 TERM2 
T1-SMT1/N 
T2-SMT2/N 
FMEAN--T1 -T2+F I 

VAR=(SMT3~NX (Tl-7 2>**2>/<N-l > 

WRITE (X, 10) 

WRITE (X, 109) 

FORMAT <11X,'T1/,7X,/T2/,7X,/ FMEAN '  , 7X , ' VAR ' ) 

WRITE  OK, 10) 

WR I TE ( X , 1 1 0 )  T1 ,T2* FMEAN , VAR 
F  ORMAT  ( 7X ,  4  ( IX , F'9 *  8 )  ) 

WRITE  OK, 10) 

STOP 

END 

FUNCTION  F  Ur.’O  <  Y ) 

REAL  FY.Y 

FUNC= ( EXP <  Y  > -1 ,0)71.718 

F 1 =FUNC 

RETURN 

ENT 

St  IP  ROUT  1 NE  RANI'KR) 

REAL  R,S,XM 
INTEGER  IS 


570“  DATA.  S,XM  / . 2397S , 823 . 53478/ 

580«  5=8#Xtf 

590s-  IS  =  S 

6  1  0=  S  =  F.r  ,  823-4617 

620-  RETURN 

630s-  END 

6  4  0  F  U  N‘  C  T 1 0  N  F  A  X  <  X  ) 

650s  F  A  >:  =  >■' 

6  70s  END 


X  X  *  X  X  *  *  X  X  t  *  X  XXX  X  *  #*#  X  x  X  *  X  X  *  XX  *  *  X  X-  X  X  X'  *  *  X  X  X;  X  *  *  *  X  *  *  *  * 
RM  TERRI  TERR 2  TERMS 

X  X  X  X  X  X  X  *  *  X  X  X  X  X  X  X  *  X  X  X  X  X  X  X  X:  X  X  X  *  X  X  X  X  X  X  X  X.  X  X  X  X  X  X  X  X  X  XXX  X 


.467170 

♦  346609 

.467170 

.014535 

<■-.  *7  r\  — 4 

4  \  . 

, OP 085 4 

,079721 

, 003465 

,630365 

.511232 

.  630365 

.014193 

.276951 

. 185740 

.276951 

.00331? 

.228132 

.  149157 

,r>pf,  •;  •7,-> 

*  0  0  £•  2  w-  / 

,024116 

,014208 

.0241 16 

.000093 

.009957 

.005825 

. 009957 

.000017 

. 349598 

.  !XC. 

,349598 

.011237 

.055349 

.033125 

.055349 

,000494 

.731006 

.626990 

,731006 

.010819 

.158201 

.099768 

.158201 

.  0034 1  4 

.  433205 

.  3 1 5595 

,433205 

.013832 

.908567 

.861913 

.908567 

.002177 

.  386255 

,274424 

,386255 

.012506 

,243985 

.160841 

.243985 

.006913 

.079153 

.047945 

.079153 

.000974 

.334626 

.231326 

. 334626 

.010671 

72559"-’ 

,620469 

,  725597 

.011052 

. 704083 

,594873 

.704083 

,011927 

.986214 

.978501 

.986214 

.00005? 

X XXX XX XXX*  XXX  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  XX  XX 
T3  T2  FMEAN  OAR 

XXXX  X X X X X X  X X X X X X X X X X X X X X X X X X X X  XXX  XXXXX  XXXXXXXXXXXX 
,356.1.494:?  .43061259  ,42553684  .00168658 
X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X  X 
61300  CM  STORAGE  USED. 

0,130  CP  SECONDS  COMPILATION  TIME. 

STOP 

15100  MAXIMUM  EXECUTION  FL . 

0.05?  CF-  SECONDS  EXECUTION  TIME. 


APPENDIX  F 


Computer  Program  Listing  with  the 
Results  for  the  Example  of  Antithetic  Variates 


L.  .  A 


1 0  0  -•  P  R  0  G  R  A  H  A  N  T  H  T  C 

iio*c  ###*#*#  *  ♦  w  *  *  .*■  *  *  *  *  *  *  #  #•  *  *  #  *  *  #  *##**)»:###»**  *  a?#***-*  *  *■  *  *  *•  *  *x  *#* 

120- C  THIS  PROGRAM  IS  AH  EXAMPLE  OF"  USING  ANTITHETIC  VARIABLES  AS 

130~ C  V ♦ R  *  T 

1  4  0  -  f ;  *  ********  #  *  *  *  -*■  *.  *  ****  *  *  *  *  *  *  *  ****  *  i :  *  *  ********  *  »  ************  * 

150-  ALPHA® «  45 

160®  WRITE < *  ,  300 > 

170*  DO  1  L«1»1C 

IDO®  8 U Ml® 0.0 

190®  SUM2®0 ♦ 0 

200®  N--20 

210®  ALP  HA® ALPHA 4  .05 

220®  HO  100  1*1, N 

230®  CALL  RANT*  < F:  1 ) 

2  AO®  7  ERM 1  *  FUNO  <  R 1  ) 

250®  F;2*l-R3 

260®  TERM2«FAX(R2) 

270®  TERM  3= ALPHA*  <  TERM! )  4  ( 1-6...PHA  )  *  (  TER  M2  ) 

2 e 0 ®  S U  Ml ®S U M 3  + T  E R M  3 

290®  100  SUM2«SUM24TERM3**'2 

300®  EI=SUM1/N 

31 0®  EI2=EI**2 

320®  S2® ( 5UM2-N*EI2  > / ( N-l ) 

330®  WRITE <  * , 200)  El, S2, ALPHA 

340®  200  FORMAT (1  OX, 'THE  ESTIMATE  *  ' , FG . 7 , 3X , ' THE  VARIANCE®' ,F8. 7, ' 

A  ®  ' 

350®  1  ,F4.2) 

360®  WRITE  OK, 300) 

370®  300  FORMAT ( 1  OX , 58 ('*')> 

380®  1  CONTINUE 

390®  END 

400®  FUNCTION  FUNC(Y) 

410®  FUNC®<E.:.XP<  Y)-l  ,0)/l  ,7182818 

420®  RETURN 

430®  END 

440®  FUNCTION  FAX(X) 

450®  FAX»X 

460®  RETURN 

470*  END 

A 8 0=  SUB R 0 U TIN E  R A N D < R ) 

490®  DATA  S,XM  /♦ 23978, 823. 53478/ 

500®  S*S*XK 

5.10®  IS®S 

520®  R®S-I S 

530®  S=R  f . 8234617 

540®  RETURN 

550®  END 


F-l 


FTN5 


61300  CM  STORAGE  USED. 

0,121  CP  SECONDS  COMPILATION  TIME. 


*  **  ******  ***** *** Hi*#:*##  ***  ******:*#*******#**************** 

THE  ESTIMATE  -  . 4627392  THE  VARIANCE*. 00042 IS  A  *  ,50 

******* * * *  * * * * * * * * * ************ * *  * ***************** * * * * * *  * 
THE  ESTIMATE  ®  .4523036  THE  VARIANCE™  .0016324  A  *  .55 

*  *  *  *  *  *  *  *  *  *  if  *  *  *  *  *  *  *  *  ******  *  *  *  *  *:  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 

THE  ESTIMATE  ®  .4281094  THE  VARIANCE--  .  0025942  A  -  .60 

*  *  *  *  *  >.  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *****  *  *  *  *;  *  *  *:  *:  *  *  *  *:  *:  *  *  *  *  *  *  *  *  *  *  *  *  *.  *:  *  * 

THE  ESTIMATE  -•  .4817804  THE  VARIANCE- . 0095206  A  «  .65 

*  *  *  *  *  *  *  *  *  *  *  *.  *  *  *  *  *  *  *  *********  *:  *  *  *.  *  *  *  *  *  *  ************  *  *  *  *  **** 

THE  ESTIMATE  ®  .4264398  THE  VARIANCE® .0111132  A  =  .70 

******  >;;  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  ********  *  *  *  *  *  *  *  *  *  *  *:  >' '  *  *  *  *****  *  *  *  *  *  *  * 
THE  ESTIMATE  ®  .4180123  THE  VARIANCE* .0257596  A  *  .75 

*  *  *  *  ***  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *.  *  *  *  *  *  *  *  a  **  *****  *  i  y ; 

THE  ESTIMATE  *  ,4686817  THE  VARIANCE® . 0356525  A  =  .80 

*  *  *  *  *  *  *  *  if  *  *  *  *  *  *  *  *  *  *  *  *  *#*  *  ******  *  *  *  *  *  *  *:  *  *.  *  *  *  *  ***  *  *  *  *  *  *  *  *  ***. 

THE  ESTIMATE  =  .4030911  THE  VARIANCE® .0403216  A  ®  .85 

* * * * *  * * * *  *  *********** * *  * * * * *  *  * * * *  * * * * * * * * *  * *  *  * * *  *  *  *  * * * *  V r ) 

■ME  ESTIMATE  «  .38718S3  THE  VARIANCE® . 0483 605  A  ®  .90 

******** m  +  *******  **************>  *  :  *******************.**** 
the:  ESTIMATE  ®  .48* 6249  THE  VARIANCE®.  073  7330  A  ®  .95 

*********  *********************************************  **** 
END  ANTHTC 

35100  MAXIMUM  EXECUTION  FT  . 

0.067  CF  SECONDS  EXECUTION  TIME. 
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